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This  project  was  concerned  with  the  development  of  biocybcrnetic  concepts 
and  methods  which  have  potential  value  for  enhancing  visual  perception  of,  and 
visual  memory  for,  scenic  materials.  The  essence  of  the  biocybcrnetic  idea  is 
that  feedback  and  control  schemes  implemented  by  various  machines  can  be  u.scd 
to  enhance  or  extend  various  aspects  of  human  performance  beyond  unaided  limits 
To  design  effective  closely-coupled  mnii-machine  systems  it  is  desirable  to  have 
the  perceptual  and  behavioral  aspects  of  human  performance  formulated  in  terms 
of  feedback  and  control  principles.  In  this  project  the  investigators  concentrijted 
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on  studying  the  feedback  and  control  t ossibilities  inherent  in  the  coupling  of  1 
visual  stimuli  to  eye-direction  and  to  the  phase  of  the  EEG  alpha  rhythm. 

They  succeeded  in  developing  state-of-the-art  systems  for  real-time  tracking  oJ 
eye-direction  and  alpha-phase;  these  computerized  tracking  systems  are  capable 
of  controlling  visual  stimuli  so  that  their  occurrence  is  conditional  upon  eye- 
direction  and  alpha-phase. 

Their  computerized  eye-tracking  system,  known  as  PERSEUS,  incorporates  the 
2-dimensional  double-Purkin je-image  eye-tracker  (DPIET)  developed  at  the 
Stanford  Research  Institute.  This  non-contacting  device  tracks  and  compares 
the  positions  of  the  first  and  fourth  Purkinje  images  formed  by  reflections  of 
infrared  light  beamed  at  the  subject's  eye.  By  comparing  the  position  of  the 
fourth  Purkinje  image  with  respect  to  the  first  Purkinje  image,  one  obtains  a 
sensitive  measure  of  eye-rotation  which  is  uncontaminated  by  eye-translation. 
PERSEUS,  which  is  implemented  in  a PDP-15/76  computer,  provides  computer- 
generated scenic  targets  for  experimental  studies,  calibration  routines  for 
the  linearization  of  eye-direction  estimates  obtained  from  the  DPIET,  and  real- 
time analysis  of  fixation  and  saccades  in  relation  to  scenic  targets.  PERSEUS 
is  capable  of  stabilizing  a computer-generated  scenic  image  upon  the  retina  of 
the  moving  eye;  it  is  also  capable  of  modifying  the  properties  computer-generated 
1 scenic  target  on  the  basis  of  current  eye-direction.  A model  for  saccadic  eye- 
movement  was  developed  which  permits  real-time  prediction  of  saccadic  destina- 
tion from  the  early  portion  of  the  saccadic  trajectory. 

'' — ^ A computerized  scheme  was  implemented  for  the  phase-contingent  analysis  of 
the  average  visual  evoked  potential. ^The  results  of  this  investigation  were 
found  to  be  compatible  with  the  hypothesis  that  cerebral  incorporation  of 
visual  sensory  inputs  is  based  upon  phase-miodulation  of  the  EEC-  alpha  rhythm. 

A more  elaborat::;  analysis  of  the  average  visual  evoked  potential  in  terms  of 
alpha  phase  probability  at  the  moment  of  photic  stimulation  also  produced 
results  consistent  with  the  phase-modulation  hypothesis.  According  to  the 
phase-modulation  theory  advanced  by  the  investigators,  the  cerebral  encoding 
of  visual  sensory  inputs  is  based  on  the  phase-difference  between  the  cortical 
alpha  phase  based  upon  autonomous  stimulation  by  the  thamalic  pacemaker  and 
the  observed  cortical  phase  which  may  be  shifted  by  the  action  of  the  sensory 
(exogenous)  input.  According  to  this  theory,  the  alpha  'carrier'  would  have 
to  be  activated  in  order  for  demodulation  of  mnemonic  playback  to  occur.  A 
nonlinear  model  of  the  human  EEG  signal  was  developed  and  tested. 

Since  practical  success  in  achieving  a biocybemetic  ' close  coupling'  of 
man  and  machine  depends  upon  an  enlighted  selection  of  suitably  efficient 
modeling  techniques,  a review  of  recursive  modeling  methods  was  undertaken  to 
provide  a systematic  classification  of  an  area  characterized  by  rapid  and 
diverse  developments. 
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FOREWORD 

The  Stanford  Biocybernetics  Project  was  begun  in  1972  with  Professor 
D.  C.  Lai  and  Professor  T.  Kailath  as  co-principal  investigators. 

Dr.  J.  Anliker  of  the  NASA/Ames  Research  Center  played  an  active  role 
as  consultant  in  the  planning  of  this  project,  made  the  facilities  ot 
his  laboratory  available  to  the  project,  and  helped  to  guide  psyclio- 
physiological  and  experimental  aspects  of  the  project.  Dr.  H.  Magnu.ski 
served  as  Post-Doctoral  Research  Fellow  during  1973-1974,  while  A.  Hua-.g, 
K.  Jacket  and  L.  Strickland  provided  programming  assistance.  Two  Ph.D. 
theses  by  J.  Nickolls  and  Arun  Shah,  both  Graduate  Student  Research 
Assistants,  were  completed  under  the  project.  After  Professor  Lai's 
return  to  Vermont  in  1975,  the  major  effort  was  contributed  by  Dr.  Anliker 
and  by  Mr.  R.  V.  Floyd  who  served  as  Scientific  Programmer.  Professor 
Martin  Morf  made  valuable  contributions  in  the  development  of  fast 
algorithms . 


T.  Kailath 
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ABSTRACT 


The  ARPA-sponsored  Biocybernetics  Project  at  Stanford 
University  has  been  concerned  with  the  development  of  bio- 
cybernetic  methods  for  the  enhancement  of  visual  perception 
of,  and  memory  for,  scenic  materials.  The  essence  of  the 
biocybernetic  concept  is  that  feedback  and  control  schemes 
implemented  by  appropriate  machines  can  be  used  to  enhance 
or  extend  various  aspects  of  human  performance  beyond  the 
unaided  limits.  In  the  present  project  the  focus  has  been 
upon  enhancing  human  perception  and  memory — in  particular, 
perception  and  memory  of  scenic  stimuli.  Since  it  is  recog- 
nized that  some  individuals  ,are  decidedly  better  visualizers 
than  others  and  since  it  is  thought  that  the  ability  to  visu- 
alize scenic  materials  would  confer  definite  advantages, 
this  project  was  conceived  to  explore  the  possibility  that 
biocybernetic  methods  might  be  employed  to  good  effect  in  the 
study  of  human  visualization  processes  and  to  determine 
whether  such  skills  might  be  trainable  through  appropriate 
feedback  or  extended  by  computerized  prosthetic  devices. 

In  this  project  we  have  concentrated  on  analyzing  the 


feedback  and  control  possibilities  inherent  in  eye-pointing 
behavior  and  in  the  electroencephalographic  alpha  rhythm. 

To  this  end  we  have  succeeded  in  developing  what  is  cuirently 
the  most  advanced  eye-tracking  system  in  existence,  and,  ir 
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collaboration  with  the  Ames  Research  Center,  we  have  also  ! 

developed  a state-of-the-art  computerized  system  for  track-  ^ 

ing  the  EEC  alpha  rhythm.  Both  of  these  systems  have  the 
capability  of  monitoring  their  biological  targets  on  an  on- 
line, real  time  basis.  In  addition,  they  can  control  visual  I 

stimulation  so  that  its  occurrence  is  conditional  upon  cur- 
rent states  in  the  biological  targets.  These  systems  are 
described  in  some  detail  in  the  body  of  this  report,  its 
appendices,  and  references.  Work  on  this  project  has  also 
been  reported  in  various  publications,  reports,  and  theses; 
two  Ph.D.  dissertations  have  been  submitted  to  the  Depart- 
ment of  Electrical  Engineering,  and  a third  Ph.D.  disserta- 
tion is  being  prepared  by  a student  in  the  Neurosciences 
Program. 

Accurate  real-time  monitoring  of  human  eye-pointing 
constitutes  a very  considerable  technical  task.  Our  advanced 
computerized  eye-tracking  system,  known  as  PERSEUS,  incor- 
porates the  2-dimensional  double-Purkinje-image  eye-tracker 
(DPIET)  described  by  Cornsweet  and  Crane  (1973)  and  updated 
by  Crane  and  Steele  (unpublished  reporc,  Stanford  Research 
Institute) . The  DPIET  is  a hardware  system  consisting  of 
an  infrared  light  source,  servo-controlled  optics,  and  asso- 
ciated electronic  circuitry.  This  device  tracks  and  com- 
pares the  positions  of  the  first  and  fourth  Purkinje  images 
formed  by  reflections  of  infrared  light  beamed  at  the  sub- 
ject's eye.  Inasmuch  as  these  two  Purkinje  images  change  | 

their  separation  when  the  eye  rotates  but  move  the  samic  amount  • 

i i i 


in  the  same  direction  when  the  eye  translates,  their  amount 
of  separation  provides  a sensitive  measure  of  eye-rotation 
which  is  uncontaminated  by  the  principal  source  of  eye-tracking 
error,  viz.,  eye-translation.  However,  the  DPIET  has  cer- 
tain limitations  which  must  be  compensated  for  by  the  compu- 
terized system.  For  instance,  the  DPIET  knows  nothing  about 
higher  order  phenomena  such  as  fixations  and  saccades;  these 
must  be  discriminated  as  temporal  patterns  appearing  in  the 
continuous  output  voltages  representing  horizontal  and  ver- 
tical eye-positions.  Nor  does  the  DPIET  have  any  provision 
for  creating  and  controlling  the  display  of  visual  stimuli. 
Finally,  there  are  the  idiosyncrasies  of  each  subject's  eye 
which  can  produce  nonlinear  distortions  in  the  eye-tracker 
estimates  of  eye-pointing  coordinates.  In  PERSEUS  the 
analysis  of  fixations  and  saccades,  the  generation  and  con- 
trol of  scenic  stimuli,  and  the  linearization  of  each  sub- 
ject's nonlinearities  are  all  handled  by  a medium-sized  disk- 
oriented  digital  computer,  the  PDP-15/76.  This  highly  ac- 
curate (less  than  4 minutes  of  arc  error  over  a field  20 
degrees  in  diameter)  eye-tracking  system  has  been  used  to 
implement  various  biocybernetic  schemes  for  controlling  scenic 
displays  on  the  basis  of  eye-position.  It  is  capable,  for 
example,  of  stabilizing  a scenic  stimulus  upon  the  subject's 
retina  without  any  attachments  to  the  eye.  While  in  this 
mode  it  can  blank  any  portion  of  the  display,  e.g.,  peri- 
phery, parafovea,  or  fovea.  PERSEUS  performs  real-time 
analyses  on  sequential  fixations  and  saccades,  numbers  scan- 


path  fixations,  and  superimposes  this  information  on  .m 


ancillary  CRT  display  of  the  visual  stimulus.  The  data  col- 
lected during  experiments  are  stored  and  subjected  to  more 
complex  forms  of  off-line  analysis.  Of  particular  interest 
are  the  visual  contact  probability  plots;  these  3-dimensional 
surfaces  are  formed  by  weighting  the  area  surrounding  each 
fixation  by  an  estimate  of  the  local  retinal  acuity.  For 
quantitative  purposes  the  contact  probability  surfaces  are 
sliced  at  equal  amplitude  intervals  to  produce  a set  of  iso- 
contours which  are  then  superimposed  upon  the  scenic  stimu- 
lus pattern.  By  comparing  these  iso-contours  with  a similar 
analysis  of  the  spatial  frequencies  in  the  stimulus  it  is 
possible  to  obtain  a measure  of  selective  attention  which 
takes  into  account  the  eye's  natural  preference  for  high 
frequency  "hot  spots." 

Our  studies  of  the  EEC  alpha  rhythm  have  been  of  two 
kinds;  (a)  detailed  off-line  analyses  of  the  individual 
sample  functions  which  are  ordinarily  subjected  to  time 
averaging  to  produce  the  visual  average  evoked  potential  and 
(b)  studies  of  phase-conditional  photic  stimulation.  For 


the  evoked  potential  studies  we  have  employed  quadrature 
analysis  for  the  definition  of  EEC  alpha  phase,  i.e.,  phase 
with  respect  to  a coherent  alpha  phase  (not  interhemispheric 
phase  differences) . In  one  study  we  examined  the  role  of 
the  contingent  alpha  phase,  i.e.,  the  phase  of  the  alpha 
rhythm  at  the  moment  of  photic  stimulation.  When  we  re- 
classified the  AEP  sample  functions  into  twenty  different 
contingent  phase  ranges,  averaged  the  sample  functions  within 
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each  phase  range  separately,  and  compared  these  contingent 
averages,  we  found  that  the  stimulus  seemed  to  phase-shift 
each  of  these  averages  in  the  latency  range  associated  with 
the  prominent  N1-P2-N2  complex  in  the  AEP , This  suggested 
that  the  stimulus  might  be  acting  to  phase-shift  the  cortical 
alpha  rhythm  away  from  the  phase  anticipated  in  the  autono- 
mous or  unstimulated  alpha  rhythm.  This  would  produce  a 
cortical  phase  at  variance  with  the  thalamic  pacemaker.  To 
check  further  into  this  possibility,  we  extracted  the  time- 
varying  phase  function  from  each  of  the  raw  EEG  sample  func- 
tions. This  gave  us  the  phase  of  the  cortical  alpha  activity 
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in  each  sample  function  at  2-msec  intervals  following  the 
stimulus.  From  these  phase  data  we  constructed  3-dimensional 
phase  probability  surfaces  (phase,  latency,  probability)  for 
each  of  the  phase-contingent  subsets  of  the  AEP  data.  Then 
we  extracted  the  modal  phase  function  from  each  of  these 
contingent  phase  probability  surfaces.  This  analysis  showed 
clearly  that  our  previous  conclusion  based  on  phase  contin- 
gent averaging  was  correct.  That  is,  the  photic  stimulus 
acts  shift  the  contingent  phases  in  the  direction  of  a "pre- 
ferred" phase  at  the  latency  corresponding  to  the  prominent 
N1-P2-N2  complex  in  the  AEP.  These  results  are  also  consis- 
tent with  our  phase  uncertainty  model  of  stimulus  encoding 
in  the  visual  cortex. 

For  our  phase-conditional  stimulation  studies  we  de- 
vised a novel  method  for  tracking  and  modeling  the  frequency 
and  phase  of  the  cortical  alpha  rhythm.  In  this  computerized 
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scheme,  the  alpha  frequency  is  monitored  continuously  by 
autocorrelating  the  EEG  signal  and  interpreting  this  trans- 
form. This  computation  is  accomplished  by  a special  purpose 
computer  (SAICOR  400-pt.  correlator) . A general  purpose  com- 
puter does  the  interpretation  into  frequency  and  uses  this 
information  to  control  the  frequency  parameter  of  a pro- 
grammable oscillator.  The  output  of  the  oscillator  (now 
tuned  to  the  alpha  frequency)  is  fed  into  another  correlator 
where  it  is  cross  correlated  with  the  current  EEG  signal. 

The  output  of  this  correlator  is  interpreted  as  phase  dif- 
ference between  the  two  signals;  this  phase  information  is 
used  to  null  the  phase  difference  by  adjusting  the  phase 
parameter  of  the  programmable  oscillator.  Through  this  com- 
plex procedure  the  output  of  the  programmable  oscillator  is 
made  to  model  closely  the  frequency  and  phase  properties  of 
the  EEG  alpha  rhythm.  From  this  point  it  is  a relatively 
simple  matter  for  the  computer  to  malce  the  presentation  of 
a photic  stimulus  conditional  upon  a specified  phase  of  the 
alpha  rhythm. 

In  a project  such  as  this,  closely-coupled  systems 
cannot  be  achieved  without  the  aid  of  efficient  estimation 
and  prediction  of  the  brain's  states  as  revealed  in  such 
outputs  as  eye-movements  and  EEG  scalp  potentials.  Although 
the  brain's  activities  are  complex  and  time-varying,  linear 
systems  modeling  techniques  can  be  used  successfully  to  pre- 
dict brain  states.  Model  parameter  estimates  can  be  updated 
in  time,  yielding  a useful  time-varying  linear  model . 
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Nonstationarities  in  observed  outputs  can  be  modeled  with 
time-varying  parameters,  or  with  a time-invariant  linear 
m.odel  with  suitable  initial  conditions.  In  the  appendix 
entitled  "A  Classification  of  Recursive  Modeling  Methods" 
we  have  presented  a systematic  classif ication  of  exact  least- 
squares  modeling  procedures  that  are  recursive  (in  model- 
order)  and  optimal  in  some  sense.  Methods  which  are  sub- 
optimal  or  approximate  are  only  briefly  indicated. 

Of  the  three  types  of  procedures  considered,  e.g., 
Riccati  or  square-root  methods,  fast  methods  utilizing  matrix 
shift-invariance  properties,  and  ladder-form  methods,  the 
latter  two  are  most  suited  to  problems  requiring  efficiency. 
In  particular,  the  recursive  (in  time)  versions  lend  them- 
selves to  on-line  or  real-time  applications  because  the 
input-output  data  are  assessed  sequentially  one  sample  at  a 
time.  Also,  new  parameter  estimates  are  available  at  each 
sample  time,  which  facilitates  the  solution  of  on-line 
control  problems. 

The  ladder-form  methods  also  have  the  desirable  prop- 
erty that  their  stability  can  be  checked  merely  by  inspec- 
tion of  the  reflection  coefficients.  In  addition,  they  are 
numerically  robust  since  the  major  operations  are  sample 
cross  correlations.  They  also  have  minimal  storage  require- 
ments for  least-squares  modeling  methods.  The  structure  of 
the  ladder-type  realization  suggests  that  the  reflection  co- 
efficients may  have  physical  significance  in  the  process 
being  modeled. 
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Finally,  it  may  be  noted  that  practical  success  in 
achieving  the  desired  close-coupling  of  man  and  machine  will 
depend  upon  an  enlightened  selection  of  suitably  efficient 


modeling  techniques.  Our  review  of  recursive  modeling 
methods  is  intended  to  provide  a systematic  classification 
of  an  area  characterized  by  rapid  and  diverse  developments. 
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INTRODUCTION 

The  purpose  of  the  Stanford  Biocybernetics  Project  was 

\ 

to  investigate  the  feasibility  of  developing  and  using 
biocybernetic  technology  (a)  to  study  visual  picture  percep- 
tion and  scenic  memory  and  (b)  to  enhance  human  memory--! n 
particular  the  quality,  persistence,  and  retrievability  of 
visual  images  of  environmental  objects.  Would  it  be  possi- 
ble, in  other  words,  to  use  real-time  com.puterized  monitors 
and  feedback  loops  to  arrange  extraordinary  coincidences  of 
eye-movement  scanpaths,  electroencephalographic  states,  and 
visual  stimulus  configurations  so  as  to  favor  the  production 
of  unusually  vivid,  persistent,  and  memorable  visual  images 
of  the  stimulus  material?  Although  the  common  adult  experi- 
ence with  post-stimulus  visual  images  is  one  of  rapid  fading, 
there  have  been  persistent  reports  that,  in  some  individuals, 
vivid  visual  images  persist  much  longer  than  is  usual  and 
that  they  can,  in  these  cases,  sometimes  be  recalled  with 
remarkable  accuracy  and  modified  at  will.  Some  investigators 
have  viewed  these  phenomena  as  evidence  of  neural  pathology 
while  others  have  tended  to  regard  them  as  manifestations  of 
unique  biological  gifts.  We  are  inclined  to  ask  whether 
they  may  be  regarded  as  indications  that  these  individuals 
simply  utilize  more  successful  strategies  for  the  impression, 
storage,  and  retrieval  of  scenic  (as  op'posed  to  verbal)  in- 
formation . 
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Available  evidence  suggests  that  unusually  strong 


i 


visual  imagery  is  most  frequently  observed  in  children  and 
only  rarely  persists  into  adulthood;  furthermore,  it  appears 
that  the  development  and  exercise  of  verbal  skills  coincides 
with  a reduced  incidence  of  exceptional  visual  imagery  of 
the  eidetic  type.  If  it  could  be  shown  that  the  basic 
capacity  for  visual  memory  is  a common  human  property  which 
is  given  a low  priority  in  a predominantly  verbal  society, 
it  might  be  possible  to  devise  suitable  training  procedures 
for  the  strengthening  and  control  of  visual  memory  as  a 
voluntary  option  in  educated  adults.  But  if,  on  the  other 
hand,  visually  gifted  individuals  should  prove  to  have 
unique  biological  endowment,  it  seems  reasonable  to  ask 
whether  computers  can  be  made  to  function  as  prosthetic  de- 
vices which  would  compensate  for  inadequate  or  missing  per- 
ceptual or  mnemonic  mechanisms,  on  the  one  hand,  or  disable 
or  suppress  normally-active  imagery-inhibiting  mechanisms, 
on  the  other  hand,  so  that  visual  imagery  can  be  made  to  be 
unusually'  vivid,  persistent,  and  retrievable.  Of  course, 
since  the  uncontrolled  perseveration  of  visual  images  is 
more  likely  to  be  distracting  rather  than  useful,  the  matter 
of  volitional  control  is  of  considerable  importance.  Another 
possibility  that  deserves  consideration  is  that,  with  compu- 
terized assistance,  even  relative  strong  visual  imagery  could 
be  further  enhanced,  thereby  extending  the  limits  of  human 
performance.  We  can  say  with  some  certainty  that  one  of  the 
principal  barriers  to  progress  in  th<?  de>siqn  of  closely-couj-iled 


PERSEUS;  A STATE-OF-THE-ART  EYE-COUPLED  SCENE  GENERATOR 

Accurate  real-time  monitoring  of  human  eye-pointing 
constitutes  a very  considerable  technical  task.  Our  advanced 
computerized  eye-tracking  system,  known  as  PERSEUS,  incor- 
porates the  2-dimensional  double-Purkin je-image  eye-tracker 
(DPIET)  described  by  Cornsweet  and  Crane  (1973)  and  updated 
by  Crane  and  Steele  (unpublished  report,  Stanford  Research 
Institute) . The  DPIET  is  a hardware  system  consisting  of 
an  infrared  light  source,  servo-controlled  optics,  and  asso- 
ciated electronic  circuitry.  This  device  tracks  and  com.- 
pares  the  positions  of  the  first  and  fourth  Purkinje  images 
formed  by  reflections  of  infrared  light  beamed  at  the  sub- 
ject's eye.  Inasmuch  as  these  two  Purkinje  images  change 
their  separation  when  the  eye  rotates  but  move  the  same 
amount  in  the  same  direction  when  the  eye  translates,  their 
amount  of  separation  provides  a sensitive  measure  of  eye- 
rotation  which  is  uncontaminated  by  the  principal  source  of 
eye-tracking  error,  viz.,  eye-translation.  However,  the 
DPIET  has  certain  limitations  which  must  be  com[>ensated  for 
by  the  computerized  system.  For  instance,  the  DPIET  knows 
nothing  about  higher  order  phenomena  such  as  fixations  and 
saccades;  these  must  be  discriminated  as  temporal  patterns 
appearing  in  the  continuous  output  voltages  representing 
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horizontal  and  vertical  eye-positions.  Nor  does  the  DPIET 
have  any  provision  for  creating  and  controlling  the  display 
of  visual  stimuli.  Finally,  there  are  the  idiosyncrasies  of 
each  subject's  eye  which  can  produce  nonlinear  distortions 
in  the  eye-tracker  estimates  of  eye-pointing  coordinates. 

In  PERSEUS  the  analysis  of  fixations  and  saccades,  the  gen- 
eration and  control  of  scenic  stimuli,  and  the  linearization 
of  each  subject's  nonlinearities  are  all  handled  by  a medium- 
sized disk-oriented  digital  computer,  the  PDP-15/76.  This 
highly  accurate  (less  than  4 minutes  of  arc  error  over  a 
field  20  degrees  in  diameter)  eye-tracking  system  has  been 
used  to  implement  various  biocybernetic  schemes  for  control- 
ling scenic  displays  on  the  basis  of  eye-position.  It  is 
capable,  for  example,  of  stabilizing  a scenic  stimulus  upon 
the  subject's  retina  without  any  attachments  to  the  eye. 

While  in  this  mode  it  can  blank  any  portion  of  the  display, 
e.g.,  periphery,  parafovea,  or  fovea.  PERSEUS  performs  real- 
time analyses  on  sequential  fixations  and  saccades,  numbers 
scanpath  fixations,  and  superim.poses  this  information  on  an 
ancillary  CRT  display  of  the  visual  stimulus.  The  data  col- 
lected during  experiments  are  stored  and  subjected  to  more 
complex  forms  of  off-line  analysis.  Of  particular  interest 
are  the  visual  contact  probability  plots;  these  3-dimensional 
surfaces  are  formed  by  weighting  the  area  surrounding  each 
fixation  by  an  estimate  of  the  local  retinal  acuity.  For 
quantitative  purposes  the  contact  probability  surfaces  are- 
sliced  at  equal  amplitude  intervals  to  produce  a set  of  iso- 
contours which  are  then  superimposed  upon  the  scenic  stimulus 


pattern.  By  comparing  these  iso-contours  with  a similar 
analysis  of  the  spatial  frequencies  in  the  stimulus  it  is 
possible  to  obtain  a measure  of  selective  attention  which 
takes  into  account  the  eye's  natural  preference  for  high 

frequency  "hot  spots."  j 

LIMITATIONS  INHERENT  IN  HARDWARE  INSTRUMENTS 

Sooner  or  later,  the  user  of  a ready-made  instrument 
must  come  to  terms  with  the  capabilities  and  the  limitations 
inherent  in  his  particular  instrument.  Efficient  exploita- 
tion of  an  instrument  requires  a good  understanding  of  the 
instrument's  input  requirements  and  its  output  limitations. 

It  is  the  user's  responsibility  to  do  the  work  necessary  for 
the  proper  application  or  the  instrument.  The  maker  of  the 
instrument  is  only  responsible  for  the  intrinsic  design  and 
construction  of  the  instrument  per  se . The  instrument-maker 
only  agrees  to  do  a limited  amount  of  the  user's  work  in 
solving  certain  kinds  of  problems.  Despite  all  of  this, 
there  is  a tendency  on  the  part  of  instrument-users  to  ex- 
pect more  help  from  the  instrument-maker  than  the  instrument- 
maker  does,  in  fact,  provide.  The  instrument-user  must  take 
up  where  the  instrument-maker  leaves  off.  He  must  learn  to 
use  the  instrument  properly  and,  when  necessary,  he  must 
undertake  the  task  of  compensating  for  the  hardware  limita- 
tions . 

The  double-Purkin je-i mage  eye-tracker  used  in  this 
project  is  no  exception  to  the  above-mentioned  rule.  The 
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The  eye-tracker  is  designed  to  operate  within  certain  limits. 
The  DPIET,  for  example,  has  a limited  angular  range  over 


which  it  can  stay  in  tracking  mode.  Basically,  this  limi- 
tation is  set  by  the  instrument's  need  to  track  the  first 
Purkinje  image  and  to  track  the  fourth  Purkinje  image  in 
relation  to  the  first.  Anything  that  interferes  with  these 
requirements  will  cause  the  instrument  either  to  lose  track 
and  so  report  or  to  lose  track  but  continue  to  report  it- 
self in  tracking  mode  (i.e.,  to  track  a spurious  image). 

It  is  the  user's  responsibility  to  arrange  conditions  which 
will  increase  the  probability  that  the  instrum.ent  will  stay 
in  the  tracking  mode  and  which  will  decrease  the  probability 
of  spurious  "tracking"  responses. 

At  the  output  end  of  the  eye-tracker,  the  instrument 
makes  no  provision  for  recognizing  when  the  eye  is  "fixating" 
or  "saccading."  Therefore,  if  this  information  is  required 
by  the  user,  he  will  have  to  do  the  work  necessary  for  ac- 
complishing this  further  classification  of  the  vertical  and 
horizontal  eye-position  data.  And,  if  his  need  for  this 
information  is  constrained  by  time  requirements,  he  may  have 
to  devise  means  of  accomplishing  this  classification  at  high 
speed.  Needless  to  say,  as  the  performance  requirements  are 
increased,  the  difficulties  proliferate. 

In  our  biocybernetic  project,  we  set  as  our  goal  the 
achievement  of  eye-coupled  control  of  scenic  displays  which 
would  permit  a scenic  tartet  to  be  "stabilized"  on  a subject's 
retina  in  a specified  location  while  allowing  the  subject  a 
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certain  amount  of  freedom  in  eye-movements.  Furthermore, 
the  scenic  content  must  be  controllable  on  the  basis  of  con- 
tingent eye-position.  That  is,  without  restraining  eye- 
movements  (within  certain  limits)  we  wanted  to  be  able  to 
modify  scenic  content  in  a specified  manner  on  the  basis  of 
current  eye-position  and,  if  possible,  on  the  basis  of 
anticipated  eye-position.  Whereas  a contact  lens-mounted 
mirror  can  achieve  image  stabilization,  the  information  con- 
cerning eye-position  is  not  made  available  for  controlling 
the  scene-generator  on  the  basis  of  eye-position;  that  is, 
since  the  mirror  moves  with  the  eye,  the  image  reflected 
off  the  mirror  will  also  move  with  the  eye  and  therefore 
image-position  on  the  retina  will  not  be  influenced  by  eye- 
movements.  Thus,  by  means  of  this  physical  attachment  of 
the  image  source  (mirror)  to  the  eye  and  the  anatomica). 
attachment  of  the  retina  to  the  eye  globe,  the  image-retina 
from  moment  to  moment  is  virtually  constant.  This  method 
does  not  compensate  for  image-retina  displacements;  instead, 
it  eliminates  their  source  by  preventing  the  slip)page  fromi 
occurring.  This  approach  eliminates  the  need  for  tracking 
and  compensating  (and  for  processing  information  essential 
for  successful  tracking  and  compensation)  but  it  requires 
an  attachment  to  the  subject's  eye,  on  the  one  hand,  and 
leaves  the  experimenter  with  the  large  and  unsolved  task  of 
assessing  eye-position.  And  without  eye-position  informa- 
tion it  is  not  possible  to  control  position-conditional 
stimulation.  Since  we  were  primarily  concerned  with  ryhcrrK-tic 
modeling  of  brain  functions,  arid  since  the  visual  system  does 
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not  solve  its  target-tracking  problems  by  means  of  mechanical 
coupling  of  the  target  to  the  eye  globe,  we  decided  to  ap- 
proach the  problem  of  visual  tracking  in  a manner  roughly 
analogous  to  that  of  the  human  visual  system.  That  is,  since 
we  know  that  the  normal  human  subject  tracks  a visual  target 
by  detecting  and  measuring  image  displacement  and  by  making 
compensatory  movements  to  move  the  target  to  a preferred 
retinal  location,  we  decided  to  play  a similar  game,  namely, 
to  monitor  changes  in  eye-position  and  to  move  the  stimulus 
in  a manner  calculated  to  control  the  placement  of  the  opti- 
cal image  on  the  retina  regardless  of  the  subject's  eye- 
movements.  If  successful  in  this  endeavor,  we  expected  to 
reap  certain  rewards  in  the  form  of  having  a unique  capacity 
to  either  assist  the  brain  in  its  tracking  efforts  or  to 
hamper  it.  On  the  one  hand,  by  assisting  the  brain  in  its 
tracking  efforts  we  might  be  able  to  extend  the  performance 
limits  characteristic  of  the  unaided  visual  system,;  on  the 
other  hand,  by  interfering  in  well-defined  ways  with  the 
subject's  tracking  responses  we  might  be  able  to  discover 
some  of  the  brain's  tracking  strategies  which  are  not  trans- 
parent to  ordinary  methods  of  experimental  analysis. 

A secondary  aim,  but  one  that  is  essential  to  develop- 
ing a sustained  assault  on  the  limits  of  biocybcrnet ic  con- 
trol, was  to  devise  a computerized  eye-movement  analyzing 

and  display-controlling  system  that  would  allow  control 
\ 

capability  to  be  continually  expanded.  By  contrast,  systen.s 
which  are  inherently  capable  of  doing  only  one  thing  ut  a 
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time  and  which  preclude  the  possibility  of  simultaneous 
processing  and  control  of  more  than  one  variable:  this  means 
that  in  order  to  mount  an  attack  on  one  front  you  must  re- 
treat on  another  front.  For  example,  if  you  had  only  the 
2-dimensional  eye-tracker  and  the  Cornsweet-Crane  optometer, 
you  have  to  decide  whether  you  wish  to  measure  eye-position 
or  visual  accommodation  because  you  could  only  use  one  in- 
strument at  a time  (although  the  optometer  has  some  capacity 
for  measuring  eye-position,  it  is  rudimentary  as  compared 
with  the  eye-tracking  capacity  of  the  eye-tracker) . The 
recently  developed  3-dimensional  eye-tracker  (Crane  and 
Steele,  unpublished  report)  overcomes  this  limitation  and 
permits  the  user  to  track  horizontal  and  vertical  eye-position 
while  simultaneously  tracking  accommodation.  As  v;e  shall 
see,  this  increase  in  the  capacity  of  the  instrument,  wliile 
it  removes  certain  experimental  limitations,  poses  many  prob- 
lems for  the  user  who  desires  to  exploit  the  possibilities 
of  the  instrument. 

THE  DOUBLE-PURKl NJE-IMAGE  EYE-TRACKER  (DPIET) 

The  two-dimensional  eye-tracker  developed  at  the  Stan- 
ford Research  Institute  (Cornsweet  and  Crane,  1973;  Crane 
and  Steele,  1976)  utilizes  an  advanced  method  of  eye-tracking 
which  is  based  on  servo-controlled  tracking  of  two  Purkinje 
images,  namely,  the  first  Purkinje  image  and  the  fourth 
Purkinje  image,  which  move  identically  during  translational 
movements  of  the  eye  but  whicdi  ntove  d i f fei  ent  i <i  1 1 y dutin<; 
rotational  movements  of  the  eye. 
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The  virtual  image  resulting  from  light  reflected  at 
the  interface  between  the  air  and  from  the  cornea  is  known 


1 

i 


as  the  first  Purkinje  image.  The  second  Purkinje  image, 
formed  by  light  reflected  at  the  interface  between  the  rear 
surface  of  the  cornea  and  the  aqueous  humor  of  the  anterior 
chamber  of  the  eye,  is  practically  coincident  with  the  first 

Purkinje  image.  The  third  Purkinje  image,  a virtual  image,  : 

is  formed  by  light  reflected  at  the  interface  between  the  | 

front  surface  of  the  lens  and  the  aqueous  humor;  this 

image  is  much  larger  and  more  diffuse  than  the  other  images  ^ 

and  is  formed,  fortunately,  in  a plane  which  is  remote  from 
the  plane  of  the  other  Purkinje  images.  The  fourth  Purkinje 
image  is  formed  by  light  reflected  at  the  interface  between 
the  rear  surface  of  the  lens  and  the  vitreous  humor  that 
fills  the  posterior  chamber  of  the  eye  globe.  The  rear  sur- 
face of  the  lens  functions  as  a concave  mirror;  its  reflected 
light  forms  a real  image  of  the  light  source. 

The  first  and  fourth  Purkinje  images  are  almost  the 
same  size,  and  are  formed  in  almost  exactly  the  same  plane. 

Due  to  the  fact  that  the  change  in  the  index  of  refraction 
at  the  lens-vitreous  interface  is  much  smaller  than  that  at 
the  air-cornea  interface,  the  intensity  of  the  fourth  Purkinje 
image  is  relatively  weak,  i.e.,  less  than  1%  of  the  intensity 
of  the  first  Purkinje  image. 

During  translatory  movement  of  the  eye  1 he  first  and 
fourth  Purkinje  images  move  in  the  same  direction  and  travel 
the  same  distance  (i.e.,  the  distance  between  them  dots  net 
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change).  However,  when  the  eye  rotates,  the  first  and 
fourth  Pur)cinje  images  move  in  opposite  directions,  causing 
them  to  move  either  closer  together  or  farther  apart.  The 
distance  separating  the  first  and  fourth  Purkinje  images 
thus  provides  a basis  for  measuring  the  angular  rotation 
of  the  eye  without  the  error  contributed  by  translatory 
movements  of  the  eye. 

The  2-dimensional  double-Purkin je-image  eye-tracker 
provides  two  separate  output  voltages  (continuous)  whose 
magnitudes  are  proportional  to  the  horizontal  and  vertical 
components  of  the  distance  separating  the  first  and  fourth 
Purkinje  images  formed  in  the  eye  in  response  to  the  input 
of  an  infrared  light  beam.  This  2-dimensional  DPIET  pro- 
vides continuous  monitoring  of  the  angular  position  of  the 
subject's  eye  with  a high  level  of  accuracy  and  without  the 
need  for  mechanical  contact  with  the  eye.  Since  the  eye- 
tracker  operates  in  the  infrared,  it  does  not  interfere  with 
normal  vision. 

The  double-Purkin je-image  eye-tracker  (DPIET)  is 
limited  in  range  by  several  factors.  First,  when  the  eye 
rotates  toward  the  infrared  input  axis  (i.e.,  to  the  right), 
the  horizontal  distance  between  the  first  and  fourth  Purkinje 
images  decreases.  When  light  from  either  the  first  or  third 
Purkinje  images  invades  the  photocell  which  is  supposed  to 
see  only  the  fourth  Purkinje  image,  it  becomes  impossible  to 
continue  to  track  the  fourth  Purkinje  image.  Second,  when 
the  eye  rotates  away  from  the  infrared  input  axis,  the 
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horizontal  distance  between  the  first  and  fourth  Purkinie 
images  increases.  Tracking  in  this  direction  becomes  impos- 
sible when  the  iris  cuts  off  the  fourth  Purkinje  image. 

Third,  vertical  rotation  of  the  eye  causes  corresponding 
vertical  separation  of  the  first  and  fourth  Purkinje  images; 
tracking  of  vertical  rotation  ceases  when  the  iris  cuts  off 
the  fourth  Purkinje  image.  In  brief,  the  range  of  the  DPIET 
is  limited  by  the  pupil  boundary  in  rotational  directions 
away  from  the  infrared  input  axis  and  by  lights  from  the 
first  and  third  Purkinje  images  confusing  the  fourth  Purkinje 

image  detector  when  they  move  too  close  to  the  fourth  Pur-  ! 

kinje  image  as  the  eye  rotates  in  the  direction  of  the  infra- 
red input  axis.  From  this  it  will  be  evident  that  pupil 

size  is  an  important  limiter  of  tracking  range;  the  larger  >; 

the  pupil  size,  the  greater  the  tracking  range  in  all  direc- 
tions except  toward  the  input  axis.  Therefore,  since  re- 
duced ambient  light  favors  a larger  pupillary  size,  it  also 
increases  the  tracking  range  of  the  DPIET. 

To  minimize  tracking  loss  during  eye  blinks,  the  gain  j 

of  each  servo  motor  driver  is  sharply  reduced  during  each  | 

I 

blink.  Blinks  are  detected  by  summing  the  outputs  from  the 

four  quadrants  of  the  first  Purkinje  photodetector.  This  i 

sum,  which  is  nom.inally  independent  of  eye  direction,  is 
constant  except  during  eye  blinks;  during  a biink  this  sirni 
will  increase  or  decrease  substantially  depending  upon  the 
amount  of  light  reflected  from  the  eyelid.  The  blink  de- 
tection circuit  indicates  that  a blink  is  in  progress  (output 
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voltage  change)  when  the  sum  is  outside  the  adjustable  upper 
and  lower  limit  values.  To  minimize  tracking  loss  during 
eye  blinks  (that  is,  to  prevent  the  tracking  servo  motors 
from  moving  arbitrarily)  the  gains  of  the  servo  motors  are 
sharply  reduced  during  each  blink.  If  the  post-blink  eye- 
position  does  not  differ  greatly  from  the  pre-blink  eye- 
position,  the  servo  controls  can  quickly  resume  tracking 
their  respective  Purkinje  images  when  the  eyelid  opens. 

Our  updated  model  of  the  DPIET  (Crane  and  Steele,  un- 
published report)  is  designed  to  tolerate  up  to  one  centi- 
meter of  variation  in  eye  position  in  horizontal,  vertical, 
and  axial  dimensions.  The  axial  tolerance  is  achieved  by 
incorporating  automatic  focusing.  The  horizontal-vertical 
tolerance  is  achieved  by  making  the  input  light  path  track 
eye  position  automatically. 

The  first  and  fourth  Purkinje  image  trackers  are  both 
provided  with  output  signals  indicating  when  they  are  un- 
locked (i.e.,  out  of  tracking  mode).  These  signals  provide 
a convenient  means  of  assessing  eye-tracking  records.  We 
have  incorporated  them  into  our  system  for  mapping  the  track- 
able  perimeter. 

PERIMETER 

In  order  to  exploit  the  full  range  of  the  eye-trackcr 
it  is  necessary  to  measure  the  limits  of  eye-tracking  in 
every  direction.  We  have  developed  a routine  for  automatic 
mapping  of  the  trackable  perimeter. 

I I 


Thi.s  routine,  called 
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PERIMETER,  allows  the  experimenter  to  select  up  to  24  radial 
lines  arranged  around  a central  field  point.  The  subject 
is  instructed  to  track  the  movement  of  the  CRT  beam  as  it 
moves  centrifugally  at  a constant  speed  (selectable  by  ex- 
perimenter) along  one  of  the  radial  lines.  The  peripheral 
tracking  limit  is  marked  when  the  fourth  Purkinje  tracker 
unlocks.  Then  the  subject  returns  to  the  central  fixation 
point  and  the  target  beam  moves  out  along  another  of  the 
radial  lines.  And  so  forth.  See  Fig.  la.  Khen  all  of  the 
selected  radii  have  been  tracked  to  tiie  limit  of  track- 
ability,  the  program  then  connects  the  adjacent  radial  limits 
to  form  a trackable  perimeter.  Each  of  the  radial  limits 
is  measured  to  the  nearest  tenth  of  a degree  and  this  value 
is  displayed  on  the  graphic  perimeter.  See  Fig.  lb. 

This  method  of  mapping  the  trackable  perimeter  allows 
the  experimenter  to  adjust  his  stimulus  materials  so  that 
they  fall  within  the  trackable  limits.  This  means  (in  the 
digital  domain)  that  a picture  can  bo  offset  vertically  and 
horizontally  (with  respect  to  the  central  fixation  point) 
until  it  is  centered  in  the  trackable  space  and  that  the 
scale  of  the  picture  can  be  increased  so  as  to  fill  the 
trackable  area  or  reduced  so  as  to  fit  inside  the  trackable 
area.  Without  accurate  measurements  of  the  trackable  peri- 
meter, the  exper im.onter  is  inclined  to  make  two  types  of 
error:  first,  in  order  to  keep  his  picture  well  inside  the 
tracking  range  of  the  eye-tracker,  he  m.ay  make  his  [)ictures 
smaller  than  necessary,  an  excessively  costly  form  of 
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insurance  against  losing  track  of  eye  movements.  Second, 
he  may  make  his  picture  so  large  that  parts  of  it  lie  out- 
side the  perimeter  and  eye  movements  directed  at  those  parts 
are  untrackable.  On  the  one  hand,  it  is  desirable  to  make 
the  picture  as  large  as  possible  in  order  to  obtain  a maxi- 
mum angular  separation  between  features;  on  the  other  hand, 
such  expansion  runs  the  risk  of  making  parts  of  the  picture 
fall  outside  the  range  of  trackability . It  is  a case  of 
the  perennial  trade-off  between  field  size  and  resolution. 
Both  the  human  eye  and  the  eye-tracking  machine  have  limited 
fields  and  limited  resolving  powers. 


LINEARIZATION  OF  EYE-TRACKER  ESTIMATES 


We  notices  that  the  response  of  the  eye-tracker  to 
the  eye-movements  of  a real  eye  was  subject  to  nonlinear 
distortions,  i.e.,  the  distortions  were  different  in  dif- 
ferent parts  of  the  trackable  space.  Mind  you,  these  non- 
linearities  are  not  large  when  compared  with  other  eye-tracking 
devices  but  they  significantly  limit  the  functional  resolving 
power  of  the  eye-tracker  for  an  extended  range  of  movements. 

The  actual  human  eye  is,  after  all,  not  the  idealized  form 
anticipated  in  the  design  of  the  hardware.  The  distortions 
observed  are  nonlinear  in  the  sense  that  no  single  set  of 
vertical  and  horizontal  gain  and  offset  adjustm.ents  will 
correct  the  errors  in  all  parts  of  trackable  sp>ace.  What  is 
required  is  local  control  of  gains  and  offsets,  i.e.,  inde- 
pendent gain  and  offset  controls  for  each  part  of  the  track- 
able  space.  By  linearization  of  these  nonlinear  distortions 


we  could  expect  to  extend  the  spatial  resolving  power  of 
our  eye-tracking  system. 

For  the  above-mentioned  purpose  we  have  implemented 
an  adaptive  filtering  scheme  to  obtain  the  required  cor- 
rections. In  our  procedure  we  ask  the  subject  to  fixate 
sequentially  each  calibration  dot  in  a matrix  of  11  x 11 
dots;  each  calibration  dot  is  located  two  degrcies  away  from, 
its  nearest  vertical  and  horizontal  neighbors.  Actually, 
the  computer  first  displays  the  entire  matrix  of  calibra- 
tion dots  (in  the  memory  mode  of  the  CRT)  and  then  brightens 
the  particular  dot  that  the  subject  is  currently  required  to 
fixate.  When  the  subject  is  satisfied  with  his  fixation  of 
that  dot,  he  presses  a button  which  causes  the  computer  to 
sample  the  X and  Y outputs  of  the  eye-tracker.  Then  the 
computer  immediately  brightens  the  next  dot  to  be  fixated, 
and  so  forth.  It  is  assigned  that  the  subject  is  able  to 
properly  fixate  the  calibration  dot  and  that  an  adequate 
selected  sample  of  the  tracker  output  will  properly  estimate 
this  fixation,  but  subject  to  the  nonlinear  distortions  we 
are  now  considering.  When  the  subject  has  completed  his 
fixation  of  all  of  the  calibration  dots,  the  computer  immedi- 
ately displays  the  matrix  of  calibration  dots  and  super  impo.stis 
the  fixation  points  as  estimated  by  the  eye-tracker.  The 
task  of  the  corrective  filter  is  to  adjust  the  eye-trackcr 
estimates  so  as  to  minimize  the  discrepancies  between  the 
known  fixation  points  (i.e-.,  the  calibration  dots)  and  the 
corresponding  eye-tracker  estimates. 
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The  developing  filters  are  graphically  displayed  in 
Figs.  2a  and  2b.  Figure  la  illustrates  the  status  of  the 
filters  prior  to  any  adaption.  The  X-f alter  surface  is  lo- 
cated in  the  lower  left  portion  of  these  figures;  the  Y- 
filter,  in  the  lower  right  portion.  The  filter  is  repre- 
sented by  a matrix  of  21  x 21  dots  {twice  the  density  of 
the  calibration  dots)  and  is  depicted  as  lying  in  a hori- 
zontal plane.  This  plane  represents  zero  correction;  dis- 
placements above  this  plane  represent  corrective  additions ; 
displacements  below  this  plane,  corrective  subtractions. 

The  upper  left  portions  of  Figs.  2a  and  2b  shov;  the  cali- 
bration points  (11  X 11) , the  eye-tracker  estimates  of  fixa- 
tions, and  ellipses  indicative  of  the  magnitude  of  the  error. 
In  the  upper  right  portions  of  these  figures  are  displayed 
three  useful  statistics  of  error;  these  are  updated  after 
each  adaptive  iteration  (an  adaptive  scan  of  all  the  data 
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samples) . The  uppermost  statistic  is  the  maximum  error  ex- 
pressed as  a percentage  of  the  viewing  angle  (in  this  case 
it  is  20  degrees) . The  second  statistic  is  the  mean  square 
error  expressed  as  a percentage  of  the  viewing  angle.  The 
third  statistic  is  similar  to  the  second  but  expressed  it. 
terms  of  CRT  (scope)  display  units  (there  being  1024  unit'; 
across  the  display).  Figure  2b  illustrates  the  condition  of 
these  displays  and  statistics  after  the  adaptive  filteiir.g 
process  has  greatly  reduced  the  nonlinear  errors  and  the 
filter  surfaces  have  approximated  their  asymptotic  values. 
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These  filters  are  somewhat  analogous  to  optical  lenses 
designed  to  correct  for  astigmatic  visual  distortions.  We 
have  observed  that  the  filter  characteristics  obtained  on 
different  days  from  the  same  subject  are  surprisingly  simi- 
lar. Stated  otherwise,  this  means  that  the  filters  are 
quasi-stationary  or  relatively  invariant  with  changes  in 
time.  Thus,  rather  than  starting  the  calibration  process 
"from  scratch"  on  each  separate  experimental  occasion,  it 
is  possible  to  install  a previously-generated  filter,  ad- 
just the  offsets  of  its  central  point,  and  adapt  from  there. 
In  the  terminology  of  control  theory,  this  usage  of  relevant 
prior  information  is  called  a feedforward  operation  because 
it  allows  the  system  to  advance  to  the  "ball  park"  of  the 
optimal  filter.  Such  exploitation  of  stationary  parameters 
of  a measurement  system  provides  for  greater  convenience 
and  greater  speed  in  the  operation  of  the  system.  In  some 
cases  exploitation  of  these  properties  is  essential  to  suc- 
cessful tracking. 

SEER:  ROUTINES  FOR  THE  CREATION,  EDITI NG , STORAGE,  AND 
DISPLAY  OF  DIGITAL  PICTURES 

In  order  to  gain  the  full  advantages  that  derive  from 
having  a real-time  monitor  of  eye-direction  it  is  desirable 
to  have  access  to  computer-generated  visual  targets.  Then 
if  the  computer  knows  where  the  eye  is  pointed,  it  can  move 
the  target  into  some  specific'd  relation  to  the  eye.  Because 
of  the  great  speed  of  the  computer,  the  visual  effect  is 


equivalent  to  physically  coupling  the  optical  image  to  the 
retina.  Yet,  unlike  stabilized  images  produced  by  optical 
systems  physically  attached  to  the  eye,  the  digital  imago 
is  subject  to  highly  controlled  manipulation.  However,  the 
generation  and  management  of  non-trivial  digital  pictures 
can  be  a rather  complex  matter. 

Conversion  from  analog  scenes  to  digital  strings  is 
facilitated  by  the  use  of  a GRAFPEN  system  which  employs 
an  ultrasound-emitting  pen  tip  and  a drawing  surface  bordered 
at  the  top  and  the  left  side  by  ultrasound  sensor  strips. 

The  GRAFPEN  system  outputs  X and  Y voltages  which  are  pro- 
portional to  horizontal  and  vertical  positions  of  the  graf- 
pen.  These  voltages  pass  to  the  A/D  converters  of  the  PDP-15 
computer.  Although  it  might  seem  that  this  arrangement  would 
solve  the  problem  of  converting  analog  scenes  into  digital 
number  strings,  it  proves  to  be  awkward  and  inefficient. 

As  our  needs  became  apparent,  we  developed  an  increasincily 
elaborate  set  of  routines  for  making  and  handlina  digital 
pictures.  We  have  used  the  acronym  SEER  (Stanford  Eye  Ex- 
periment Routines)  to  designate  the  software  developed  fci' 
this  purpose. 

The  work  of  making  and  using  digital  [.'icturc-s  can  be 
divided  conveniently  into  three  parts.  First,  there  is  t )ie 
work  of  the  copyist  who  has  to  trace  the  analog  picture  con- 
tours with  the  grafpen  and  control  the  com.puter's  samj  ling, 
labelling,  and  storing  operations.  St-cor-.d,  theie  is  the 
task  of  assembling  contour  sinjment s into  a "picture,"  editing 


it,  naming  it,  and  storing  it  as  an  entity.  Third,  there  is 
the  task  of  assembling  a sequence  of  pictures  for  an  expert  - 
ment.  The  various  SEER  routines  arc  designed  to  facilitate 
this  work.  The  transfer  of  these  digital  pictures  from  one 
type  of  storage  to  another  must  also  be  provided;  that  is, 
the  digital  points  have  to  be  transferred  from  core  memory 
to  fixed-head  disk,  or  from  fixed-head  disk  to  DECtape  or 
9-track  digital  tape,  or  from  fixed-head  disk  to  movable- 
head  disk,  or  in  the  reverse  direction.  The  various  SEEP 
routines  are  designed  to  make  these  transfers  relatively 
easy . 

In  most  of  these  activities  it  is  desirable  that  the 
experimenter  be  able  to  see  what  he  is  working  with.  There- 
fore, we  have  provided  extensive  display  capabilities. 

Now,  to  consider  the  work  of  the  copyist:  the  scenic 
material  to  be  copied  is  rear-projected  onto  the  ground  glass 
surface  of  the  grafpen  drawing  area.  The  copyist  must  then 
use  his  eyes  to  guide  his  hand  as  it  traces  a selected  con- 
tour with  the  tip  of  the  grafpen.  Since  this  form  of  copy- 
ing is  tedious  at  best,  it  is  convenient  to  arrange  the  sys- 
tem so  that  he  has  flexible  command  of  the  computer.  The' 
copyist's  dialogue  with  the  computer  is  carried  out  through 
the  use  of  the  grafpen.  A computer  command  "menu"  is  pro- 
vided in  the  lower  right-hand  corner  of  the  grafpen  drawing 
surface.  By  touching  the  grafpen  to  the  various  designated 
areas,  the  copyist  can  request  that  the  com.puter  give  him  an 
element  number;  the  computer's  selected  number  is  then  displayed 
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by  the  computer  on  a display  CRT  located  near  the  drawing 
surface.  In  the  MENU-mode  the  copyist  can  draw,  erase, 
label,  and  store  up  to  1000  picture  segments  at  a tiire. 

We  shall  not  elaborate  on  the  available  menu  commands 
other  than  to  mention  that  in  addition  to  the  more  obvious 
computer  commands  we  have  provided  some  "geometrical  assis- 
tance" commands  which  allow  the  copyist  to  enlist  com.puter 
assistance  in  drawing  straight  lines,  circles,  arcs,  and 
rectangles.  Also,  to  economize  on  digital  storage  and  to 
obtain  an  approximately  constant  density  of  sample  points 
from  a picture  contour,  there  is  a spatial  sampling  routine 
which  stores  a new  data-point  only  if  it  is  located  seme 
specified  distance  away  from  the  previously-stored  data- 
point. 

Equalizing  the  sample-density  of  the  picture  contours 
not  only  economizes  on  storage  but  it  lays  the  foundation 
for  meaningful  manipulation  of  contour  density.  Thus,  for 
example,  it  is  possible  to  get  reasonably  good  tracings  by 
displaying  every  other  point,  or  some  suitable  fixed  ratio 
of  points.  Or,  if  you  wish  to  generate  a moving  target 
point  for  the  eye  to  follow,  a fixed  rate  of  display  of  the 
points  will  appear  to  speec  up  in  approaching  curves  and 
corners  and  will  appear  to  slow  down  in  the  straight  sec- 
tions. By  increasing  the  sample  density  in  curves  and 
corners  it  is  possible  to  equalize  the  apparent  speed  of 
the  moving  light  beam  on  the  display  CRT.  Since  we  could 
not  arrive  at  any  generally  satisfactory  formula  that  would 


produce  subjective  equality  of  trace-speed,  we  hit  upon  an 
empirical  solution  whereby  the  editor  is  permitted  to  select 
any  particular  segment  by  marking  two  limiting  points  on  the 
contour;  then  he  can  draw  with  the  grafpen  whatever  function 
seems  reasonable  and  the  density-adjusting  routine  will  in- 
crease or  decrease  the  sample  density  along  that  segment  in 
accordance  with  the  selected  function.  Then  the  editor  chn 
request  that  the  segment  be  traced  on  the  CRT  using  the  ad- 
justed density.  If  he  is  dissatisfied  with  the  result,  he 
can  modify  the  density  function  and  see  if  that  produces  a 
better  subjective  effect;  if  he  is  satisfied  with  the  re- 
sult, he  can  move  on  to  the  next  segment.  Finally,  we  should 
mention  that  the  equal-density  contour  drawings  provide  the 
basis  for  generating  the  spatial  density  matrices  for  these 
pictures  (to  be  described  later) . 

Having  collected  a set  of  picture  segments,  the  copyist 
can  then  proceed  with  the  task  of  assembling  "pictures"  from 
picture-segments.  He  can  edit  his  picture  in  many  ways. 

For  instance,  he  can  magnify  it  or  minify  it,  he  can  offset 
it  on  horizontal  and  vertical  axes,  and  he  can  rotate  it  (yaw, 
pitch,  roll) . Then  he  can  give  this  edited  version  of  the 
picture  a label  before  storing  it  in  one  of  several  convenient 
places . 

Finally,  he  can  assemble  a set  of  "pictures"  into  an 
experimental  sequence,  specifying  various  display  param.ctei  s 
for  the  experiment.  With  t:hese  manipulative  capabilities 
it  is  possible  to  utilize  digital  pictures  much  the  way  that 


one  would  use  a set  of  photographic  slides  and  a projector. 
However,  being  in  the  digital  domain,  the  digital  picture 
can  be  manipulated  by  the  computer  in  some  extraordinary 
ways.  On  the  negative  side,  there  is  the  restriction  (in 
our  set-up,  at  least)  to  monochromatic  contour-pictures  or 
"cartoons . " 

As  a future  development,  we  foresee  the  coordination 
of  visually-rich  color  T.V.-type  of  display  with  the  sche- 
matic contour-pictures,  allowing  photographs  to  be  overlaid 
by  corresponding  schematics  on  a frame-by-frame  basis. 

FIXATION-CONDITIONAL  STIMULUS  PRESENTATION 

In  many  experiments,  the  subject  is  asked  to  fixate  a 
centrally  located  target  as  a pre-condition  for  the  exposure 
of  the  stimulus  pattern.  From  the  point  of  view  of  experi- 
mental design,  fixation  of  the  central  point  is  a stimulus 
parameter  in  the  sense  that  the  specified  fixation  arranges 
for  consistent  retinal  placement  of  the  stimulus  pattern. 

Yet  it  is  seldom  that  any  precautions  are  taken  to  monitor 
or  control  this  parameter  other  than  to  instruct  the  subject 
to  fixate  and  expect  his  cooperation.  A more  rigorous  ap- 
proach would  be  to  arrange  for  electronic  monitoring  of  eye- 
pointing and  for  electronic  circuitry  for  making  stimulus 
oresentation  conditional  upon  fixating  within  prescribed 
spatial  boundaries  for  some  prescribed  period  of  time.  This 
arrangement  has  the  advantage  that  the  stimulus  exposure  can 
be  extended  as  long  as  the  subject  keeps  his  eye  within  the 


prescribed  spatial  boundaries;  the  alternative,  of  course, 
is  to  employ  tachistoscopic  exposures  which  are  sufficiently 
brief  as  to  preclude  significant  reactive  eye-movements.  A 
computerized  fixation-conditional  stimulus  controller  is 
illustrated  in  Figs.  3 and  4.  This  is  actually  the  simplest 
use  of  a more  general  system  of  visua]  scanpath  analysis 
which  is  described  below. 

REAL-TIME  SCANPATH  ANALYSIS 

Automatic  methods  for  the  off-line  analysis  of  scan- 
paths  are  described  in  Shah  (1977)  . ^^e  shall  now  describe 

a real-time  method  for  the  analysis  of  scanpaths.  It  is  evi- 
dent that  an  effective  real-time  scanpath  analyzer  is  an 
essential  development  in  the  process  of  implementing  cyber- 
netic models  for  the  prediction  and  control  of  ongoing  eye- 
pointing behavior. 

This  computer  program  was  first  implemented  at  Ames 
Research  Center  on  an  SEL-840  computer  which  iias  an  flvans 
and  Sutherland  display  system.  Two  large  CRTs  are  utilized: 
(1)  a monitor  scope  on  which  are  displayed  the  various 
analyses  that  are  of  interest  to  the  experimenter,  and  (2) 
a stimulus  display  CRT  which  is  viewed  by  the  experimental 
subject.  During  the  past  year  we  have  dev'eloped  an  equi- 
valent real-time  scanpath  program  for  our  PDP-15/7fi  at  Stan- 
ford . 

The  subject  is  seated  in  front  of  the  stimulus  CRT, 
his  head  stabilized  by  a bitebar,  with  the  distance  Letwcer 
the  eye  and  the  CRT  surface  adjusted  so  that  12  inches  in 
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either  the  vertical  or.  horizontal  dimension  on  the  scope 
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corresponds  to  20°  of  visual  angle.  His  right  eye  is  moni- 
tored by  a 2-dimensional  eye-tracker  which  provides  output 
voltages  proportional  to  horizontal  and  vertical  eye- 
direction  . 

Following  completion  of  suitable  calibration  proce- 
dures, the  subject  is  presented  with  a cartoon-type  stimu- 
lus (out  of  digital  storage)  on  his  display  scope  and  in- 
structed to  examine  it  freely  for  a specified  number  of 
seconds  or  for  a specified  number  of  fixations,  after  which 
the  cartoon  disappears  and  the  central  fixation  spot  is 
restored.  While  the  subject  is  examining  the  cartoon,  the^ 

X and  Y voltages  from  the  eye  monitor  are  entered  into  X and 
Y voltage  distributions  located  below  and  to  the  left,  res- 
pectively, of  the  fixation  circle  (see  Fig.  3).  The  entries 
in  these  distributions  are  cumulated  in  accordance  with  a 
sliding  time-window  (the  width  of  the  window  being  speci- 
fied in  milliseconds  by  the  experimenter) . The  peak  of  each 
distribution  (X  and  Y)  is  treated  as  the  best  estimate  of 
the  eye-pointing  during  the  time  constant  characterizing 
the  width  of  the  sliding  window.  This  X,Y  location  is  made 
by  the  center  of  a circle  (solid  line)  the  radius  of  which 
is  controllable  by  the  exper im.enter . We  sometimes  refer  to 
this  circle  as  a "fixation  circle"  or  as  a "target  circle" 
depending  upon  the  experimental  use  being  made  of  it.  1 f v'e 
set  minimum  amplitude  criteria  for  the  distribution  peaks, 
we  have  a basis  for  monitoring  the  "swell  tinie”  or  "spatial 
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invariance"  as  a function  of  time,  this  is  a reasonably 
eff€?ctive  way  of  defining  the  beginning  of  each  fixation. 

The  method  for  defining  the  end  of  a fixation  is  described 
in  the  next  paragraph. 

A second  circle,  sharing  the  same  center  as  the  tar- 
get circle,  has  an  independently  variable  radius  which  is 
equal  to  or  greater  than  that  of  the  target  circle.  This 
circle,  which  appears  on  the  monitor  scope  as  a dotted  line 
(see  Fig.  3) , is  used  to  set  a spatial  variance  limit  around 
the  target  circle.  When  the  instantaneous  eye-position  ex- 
ceeds the  perimeter  of  the  variance  circle,  it  is  assumed 
that  the  eye  is  making  a saccadic  movement,  or  at  least,  that 
it  is  moving  sufficiently  far  away  that  the  information 
accumulated  in  the  distributions  is  irrelevant,  so  the  dis- 
tributions are  re-set  to  zero  and  the  existing  target  circle 
vanishes.  When  a new  fixation  is  detected,  a new  target 
circle  and  a new  variance  circle  will  appear.  And  so  forth. 

Another  important  feature  of  this  program  is  the  arrange- 
ment whereby  if  the  instantaneous  eye  position  (the  small 
square  inside  the  target  circle  in  Fig.  3)  remains  inside 
the  target  circle  for  a specified  nu nber  of  milliseconds 
(determined  by  the  experimenter),  a voltage  appears  on  an 
appropriate  analog  output  line  and  the  word  ON  is  displayed 
on  the  monitor  scope  until  the  instantaneous  eye  position 
travels  outside  the  target  circle;  v;hen  the  output  v'oltage 
goes  to  ground  the  word  OFF  replaces  the  word  ON  (see  Fia. 
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This  ON/OFF  feature  provides  the  basis  for  the  fixation- 
conditional  visual  stimulation  described  earlier.  The  real- 
time scanpath  program  provides  for  the  collection  of  up  to 
100  fixation-locations  before  it  must  pause  momentarily  to 
transfer  these  data  to  another  storage  area.  In  the  scan- 
path  display,  a cross  appears  at  each  fixation-location  and 
each  cross  is  numbered  in  accordance  with  its  ordinal  posi- 
tion in  the  scanpath. 

Figure  5a  and  Fig.  5b  contain  a cartoon  of  a paint- 
ing which  depicts  Judith  with  the  head  of  Holofernes.  In 
Fig.  5a  the  current  location  of  the  point-of -regard  is  marked 
by  the  small  circle  at  the  tip  of  Holofernes'  nose.  This 
eye-position  is  indicated  by  the  X-probability  and  Y-probabi 1 ity 
distributions  located  below  and  to  the  left,  respectively. 

When  these  probability  distributions  reach  some  selected  peak 
amplitude,  the  computer  recognizes  that  location  as  a fixa- 
tion and  deposits  a small  cross  in  the  display  to  mark  that 
spot.  Of  course,  information  concerning  the  location,  dura- 
tion, and  scanpath  number  is  stored  for  subsequent  use. 

We  have  employed  a nun±>er  of  different  techniejues  for 
the  analysis  of  scanpath  data.  One  of  the  current  develop- 
ments is  illustrated  in  Figs.  6a  and  6b.  Although  the  usual 
endpoint  of  fixation-detections  is  reached  when  fixation 
locations  are  superimposed  on  the  scenic  target,  it  is  j^ossi- 
ble  to  consider  more  complex  transforms.  In  this  appioach, 
we  have  tried  to  give  recognition  to  the  fact  that  some  inf  t- 
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FIGURE  7 

The  Iso-contours  of  Fi^.  6b  suporimposcd  on  the  pictui-c  of  Juclilh  with  t ho 
ljuuci  oj'  llolofornos.  Obviously,  Lho  spread  ol  visual  oontai-l.  around  llio 
contors  of  attontlon  wifi  dopund  ufjon  Lho  shaiKs  of  lho  oontact  probability 
woiKhts  ontcred  around  each  fixation  point.  This  system  is  flexible  in  that 
:L  allows  Lho  expo riiiion to r to  seloct  his  weif^htiiif;  pattern  on  whatever  basis 
ho  <looms  advisable.  I’rosuinably , the  pattern  of  eontaot  probability  ehani;es 
depending  upon  envi ronmonlal  illumination,  subject  motivation,  etc.  We  would 
suppose  that  tho  retinal  acuity  function  ]jrovi<les  good  basic  weighting  which 
could  be  modified  on  tho  basis  of  other  ciinside i-a t ions . 
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through  the  foveal  portion.  In  this  case  wc  have  made  a 
weighted  entry  into  every  portion  of  the  picture  space  sur- 
rounding the  fixation  point,  based  upon  some  acuity  estimate. 
The  exact  form  of  the  acuity  weighting  is  imjn.ater ia  1 for  the 
current  exercise;  the  program  permits  the  entry  of  any  speci- 
fied set  of  weights.  What  this  form  of  display  helps  to 
formulate  is  some  concept  which  we  might  call  "visual  con- 
tact probability."  Figure  6b  illustrates  a sample  probability 
matrix  (simulated  3-dimensional  plot).  Figure  6a  shows 
(enlarged)  the  equal  amplitude  transform  (or  "geophysical 
map")  of  the  data  in  Fig.  6b.  Figure  7 illustrates  such  an 
equal  amplitude  transform  of  the  contact  probability  matrix 
for  a set  of  fixation  obtained  from  a subject  who  was  given 
an  opportunity  to  examine  this  scene.  The  equal  amplitude 
contours  are  projected  onto  the  cartoon  (Judith  with  the 
head  of  Holofernes) . /Note:  the  lower  right  set  of  contours-- 
off  the  picture--resulted  from  the  fact  that  the  subject  left 
the  bitebar  before  sampling  was  terminated;  that  "fixation" 
is  the  default  location  of  the  tracking  system  in  the  abs€;nce 
of  a trackable  eye./  A dynamic  version  of  this  contact 
probability  display  provides  a view  of  the  contact  probability 
for  a relatively  narrow  sliding  window  (of  tim.o)  but  this 
display  is  more  suited  to  motion  picture  representation. 

MODELING  SACCADIC  EYE -MOVEMENTS 

A number  of  models  have  been  proposed  in  the  litera- 
ture for  the  saccadic  control  system.  While  there  is  a 


general  agreement  regarding  the  various  characteristics 
discussed  above,  there  is  no  single  model  which  explains 
them  all. 

A natural  division  of  the  system  is  to  have  a controller 
driving  a plant.  The  controller  is  the  brain  with  the  target 
signal  as  input.  The  nerve  signals  produced  by  this  control le 
drive  the  plant  which  is  composed  of  the  eye  muscle  assembly. 
Models  have  been  proposed  for  both  parts.  Westheimer  (1959) 
suggested  a second-order  linear  model  for  the  plant.  When 
such  a system  is  driven  by  a step  input,  a saccade-like  tra- 
jectory results.  In  Westheimer 's  model,  the  plant  is  under- 
damped. The  model  parameters  have  no  physiological  signifi- 
cance. Yarbus  (1967)  suggested  that  the  saccadic  velocity 
curve  is  sinusoidal.  This,  however,  is  an  over-simplification 
Based  on  the  200  ms  reaction  time  to  a pulse  target  movement. 
Young  and  Stark  (1962)  proposed  a sampled  data  control  system 
which  samples  the  error  and  then  executes  a correcting  move- 
ment 200  ms  later.  The  target  movement  in  this  interval 
does  not  cause  this  initial  response  to  change. 

The  sampled  data  model  was  presented  by  Young  and 
Stark  (1962)  and  analyzed  in  detail  by  Young  (1962).  It  is 
an  elegant  app’lication  of  classical  control  theory  to  the 
study  of  a physiological  system..  For  the  plant,  this  model 
used  Westheimer 's  (1959)  underdamped  second-order  system. 

The  model  is  adequate  for  simple  saccadic  responses  tut  can- 
not predict  the  correct  answer  to  a pulse-stei  tnrc;et.  Tl'.e 
controller  has  to  be  more  complex  to  product-  the  variet-^-  ■ *■ 


responses  observed  for  different  inputs.  Actually,  as 
Robinson  (1973)  pointed  out,  the  plant  is  overdamped  and 
needs  a pulse  superimposed  on  a step  (a  pulse-step)  as 
input  to  produce  the  observed  saccadic  trajectory.  Robin- 
son found  evidence  that  the  force  applied  to  the  extraocular 
muscles  during  a saccade  is  actually  a pulse-step.  He  used 
data  gained  through  external  loading  on  the  eyeball  as  well 
as  through  experiments  under  normal  conditions  to  determine 
the  parameters  of  the  second-order  plant  m.odel.  Robinson 
(1973)  has  suggested  a series  of  models  all  using  the  sam.e 
plant  but  with  controllers  of  increasing  com.plexity. 

Some  of  the  shortcomings  of  Robinson's  models  have 
been  pointed  out  by  Cook  and  Stark  (1967) . To  understand 
these,  we  must  refer  to  some  neurophysiological  conceptions 
of  the  saccadic  system.  Very  briefly,  the  signals  from  the 
retina  reach  area  17  of  the  occipital  lobe  via  the  optic 
nerve  and  the  lateral  geniculate  body.  The  error  signal  is 
probably  formed  in  areas  18  and  19  and  processed  to  form  the 
saccadic  motors  command  signal  in  the  frontal  eye  fields 
(Fuchs,  1971).  This  command  signal  ultimately  reaches  the 
oculomotor  nuclei.  The  force  applied  to  the  eyeball  is  a 
result  of  the  firing  of  the  oculomotor  nuclei  and  the  action 
of  the  extraocular  nuscles.  Robinson  pointed  out  that  the 
firing  pattern  of  these  nuclei  follows  a pulse-step,  with  a 
high  firing  rate  during  the  saccade  which  settles  down  to 
a steady  value  after  the  saccade.  To  produce  this  pulse-step, 
he  proposed  that  the  error  signal  drives  a pulse  generator 


which  produces  a pulse,  the  area  under  which  equals  rhe  de- 
sired saccade  amplitude.  This  pulse  is  fed  to  a network 
consisting  of  a neural  integrator  in  parallel  with  a feed- 
forward path  (the  medial  longitudinal  fasciculus).  The 
output  of  this  network  represents  the  firing  rate  of  t)ie 
oculomotor  neurons.  It  then  drives  the  second-order  over- 
damped  plant  (i.e,,  the  extraocular  muscles  and  attached 
eyeball)  to  produce  the  saccadic  eye-movement. 

There  is  considerable  physiological  evidence  indicating 
the  presence  of  a neural  integrator.  Robinson  (1972)  has 
combined  the  vestibul o-ocular  system,  the  smooth  pursuit 
system,  the  saccadic  system,  and  the  vergence  sy.stem.  into  a 
single  model  sharing  a neural  integrator  and  a feed-forward 
path.  Also  as  pointed  out  above,  the  plant  parameters  used 
in  his  model  were  derived  from  physiological  experiments. 
Thus,  it  is  clear  that  Robinson's  model  and  its  parameters 
have  considerable  physiological  significance,  which  is 
essential  in  a model  that  purports  to  explain  a biological 
system. 

Robinson's  model  does  have  some  limitations,  however. 
The  eyeball  is  controlled  by  3 pairs  of  extraocular  m.uscles 
used  for  rotation  about  the  X,  Y,  and  Z axes  respectively. 

For  example,  during  a horizontal  saccade,  one  mus.cle  (thc' 
agonist)  shortens  in  response  to  an  increase  in  neural 
activity  whereas  the  other  muscle  (the  antagonist)  lengthens 
in  response  to  co-respond j ng  decreasf-  in  neural  activity. 

Cook  and  Stark  (1967)  have  pointed  out  th..t  Rc'bir.son  does 
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not  describe  what  portion  of  the  driving  force  is  attributed 
to  the  agonist  and  what  portion  to  the  antagonist.  The 
model  is  incomplete  in  that  sense.  They  also  reported  that 
the  velocity  curves  of  Robinson's  model  do  not  agree  with 
experimental  data. 

By  considering  muscle  dynamics  and  the  different  forces 
applied  to  the  agonist  and  the  antagonist  muscles.  Cook  and 
Stark  (1967)  proposed  a nonlinear  model  for  the  plant.  The 
model  parameters  were  determined  from  physiological  experi- 
ments. They  also  published  comparisons  between  (a)  model- 
simulated  eye-position  and  velocity  curves  and  (h)  experi- 
mental data  (Cook  and  Stark,  1968).  Clark  and  Stark  (1976) 
reported  that  this  model  firoduces  time-optimal  responses. 

In  its  full  form,  the  model  proposed  Ly  Cook  and  St  irk 
is  a sixth-order  nonlinear  system  and  is  quite  comj  lex:  tiius 
it  is  clear  that  the  saccadic  system  cannot  be  fully  n t r-  - 
sented  by  any  simple  model.  The  separation  into  a rtr  . r 
and  a plant  is  itself  a simplification  which  may  r.'t  1>  t u ■. 
correct.  For  full  understanding  of  the  :..n-'-.dio  ',st>r,  ‘ 
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and  prediction.  Because  Robinson's  model  is  founded  on 
considerable  neurophysiological  evidence  and  because  of  its 
relative  simplicity,  we  have  used  it  as  a starting  point 
for  our  model  development.  We  have  improved  upon  his  model 
so  that  we  obtain  better  fit  of  the  model  responses  with 
experimental  data.  Our  model  parameters  are  estimated  by 
curve  fitting  techniques.  As  a result,  we  have  developed 
a model  that,  although  it  may  not  be  completely  sound 
physiologically,  does  in  fact  reproduce  the  input-output 
behavior  of  the  saccadic  system  with  sufficient  accuracy  to 
be  adequate  for  the  purpose  of  monitoring  and  prediction. 

Let  us  look  at  some  of  the  shortcomings  of  Robinson's 
model  and  discuss  our  modification  for  remedying  the  situ- 
ation. Firstly,  it  is  quite  clear  that  the  saccadic  system 
is  a nonlinear  system,  as  described  before;  yet  we  are  usinc 
a linear  model.  This  use  is  justified  partly  by  the  fact 
that  the  controller,  which  generates  the  input  to  tlie  plar.t 
in  response  to  the  visual  stimulus,  is  nonlinear,  and  partly 
by  the  need  for  computational  simplicity.  The  magnitude  of 
the  visual  stimulus  affects  both  the  height  .ind  width  of 
the  pulse  input  to  the  plant  so  that  even  though  the  model 
itstjlf  is  linear,  velocity  saturation  and  depend<*nce  of 
saccade  duration  on  amplitude  arc  observed.  Within  the 
plant  itself,  the  nonlinear  m.uscle  dynamics  have  been  apf.'roxi- 
mated  by  a linear  system. 

Cook  and  Stark  (1967)  pointed  out  that  the  velocity 
curves  produced  by  Robinson's  model  did  not  m.itch  experimentally- 
observed  velocity  curves.  In  addition,  the  model  as  initially 


presented  does  not  produce  overshoots  or  undershoots  in  the 
saccade.  Robinson  (1973)  suggested  a slight  change  to  one- 
of  the  time  constants  to  produce  the  desired  response.  How- 
ever, we  found  that  such  a modification  still  does  net  match 
the  observed  overshoots.  This  is  a serious  defect,  particu- 
larly since  we  are  interested  in  the  input-output  character- 
istics and  not  the  physiology  per  se . We  have  made  an  impor- 
tant m.odif ication  to  Robinson's  model  by  changing  the  form 
of  the  input  to  the  plant.  Interestingly  enough,  this  is 
the  form  of  the  input  that  would  result  from  Cook  and  Stark's 
model  if  we  ignored  the  differences  between  the  agonist  and 
the  antagonist  muscles  and  simply  added  (algebraically)  the 
forces  on  the  two.  This  change  enables  us  to  accurately  fit 
the  model  generated  position  and  velocity  responses  to  the 
experimentally-observed  ones  for  different  types  of  saccades. 
This  improves  upon  Robinson's  model  by  removing  its  main 
distraction . 

Until  now,  most  models  of  the  saccadic  system  have 
attempted  to  explain  only  the  gross  behavior  of  the  system. 
Essentially,  that  means  a study  of  the  static  behavior  of 
the  system,  with  the  primary  emphasis  being  on  whether  or 
not  a saccade  will  be  produced  in  response  to  a specific 
stimulus.  The  detailed  dynamic  characteristics  of  the  sac- 
cade have  been  of  secondary  importance.  Only  Cook  and  Stark 
have  attempted  to  fit  model-generated  responses  to  actual 
ones.  We  have  essentially  used  an  engineering  approach  to 
study  saccadic  movements.  By  using  the  metficids  of  parameter 
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estimation  and  prediction  we  have  produced  a model  which  is 
remarkably  accurate  in  its  tracking  and  prediction  of  sac- 
cadic eye  movements.  A detailed  account  of  this  model  is 
given  in  a Ph.D.  dissertation  by  Arun  Shah  (1977),  a copy  of 
which  is  submitted  with  this  report. 
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ALPHA  PHASE  CONTINGENCY  ANALYSIS  OF  THE  AVERAGE  VISUAL  EVOKED 


POTENTIAL 

Although  phase  analysis  of  the  KEG  alpha  rhythm  has  i 

interested  a number  of  neuroscientists,  there  aie,  to  our  | 

i 

knowledge,  no  published  accounts  of  successful  system, atic  i 

analyses  which  relate  alpha  phase  contingency  to  the  tradi-  | 

tional  average  visual  evoked  potential  (AVEP) . One  approach  | 

to  the  analysis  of  alpha  phase  contingency  effects  which  has  I 

been  adopted  by  a number  of  investigators  (Walter  ^ 1946;  j 

Walter  and  Walter,  1949;  Turton,  1952;  Bekkering  and  Storm  van 
Leeuwen,  1954;  Bechtereva  and  Usov,  1960;  Callaway  and  Yeager, 

1960;  Callaway,  1961;  Bechtereva  and  Zontov,  1 962  ; Dustm.an  and  ] 

Beck, 1965)  has  been  to  use  what  we  shall  call  "phase  condi-  , 

1 

tional  stimulation  techniques"  which  simplify  problems  in 

the  collection  and  management  of  data  by  arranging  to  stimu-  j 

late  only  at  specified  alpha  phases.  This  approach  )ias  the 
advantage  of  requiring  relatively  modest  instrumentation 
inasmuch  as  only  one  class  of  phase  contingent  sam.ple-functions 
is  collected  and  averaged  at  a time.  Intermediate  approaches 
which  provide  for  the  analysis  of  more  than  one  phase  con- 
tingency class  have  been  developed  by  Callaway  and  Layne  i 

(1964)  who  arranged  to  stimulate  and  collect  data  for  10  dif- 
ferent ranges  of  alpha  p'hnse,  and  fiy  Per'ond  and  Lesevre  (l‘^6I) 
who  investigated  4 different  i anues  of  ilpha  phasi'.  The 
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disadvantage  of  these  phase  conditional  stiniulation  tech- 
niques stems  primarily  from  the  limitations  imposed  by  the 
demand  for  real-time  monitoring  of  the  phase  variable  and 
for  real-time  decisions  about  delivering  or  withholding 
stimuli.  But  possibly  of  greater  importance  is  the  fac* 
that  the  abandonment  of  the  unconditional  or  noncontingent 
pattern  of  stimulation  in  favor  of  the  phase  conditional 
pattern  also  resulted  in  the  production  of  data  which  are  not 
easily  and  directly  compared  with  traditional  AVER  data. 

In  contrast  with  these  approaches,  we  (Anliker  and 
Floyd,  1977a)  have  elected  to  expand  the  traditional  AVER 
analysis  by  implementing  a completely  off-line  phase  contin- 
gency analysis  of  the  sample-functions  from  which  an  AV''EF 
can  be  extracted.  In  this  way,  we  gain  certain  fairly  ob- 
vious analytical  advantages  from  the  use  of  relatively  power- 
ful, but  time-consuming,  digital  filtering  and  classification 
algorithms.  And  we  obtain  results  which  can  be  directly 
compared  with  the  AVER.  Nevertheless,  we  should  point  out 
that  our  approach  has  one  major  drawback,  viz.,  it  entails 
the  acceptance  of  a relatively  large  computational  buider. 

EXPERIMENTAL  RROCEDURE 

The  EEC  data  were  collected  from  ? healthy  normal  young 
men  who  exhibited  prominent  alpha  rhythm.s.  l'lectrf)de  coup- 
lings were  "monopolar,"  the  active  electrode  beinq  placed 
over  the  occipital  region  of  the  scalp  ( T err  from  the  midi  i nr 
and  2 cm  above  the  inion)  and  the  inact  ive  icftiencf  di-rived 


from  yoked  electrodes  clipped  to  the  right  and  left  ear 
lobes.  The  ground  electrode  was  placed  over  the  mastoid 
process.  Right  and  left  occipital  responses  were  recorded 
simultaneously  through  amplifiers  in  a Grass  model  78  poly- 
graph onto  separate  channels  of  a Honeywell  7600  instrumem- 
tation  tape  recorder.  A Grass  PS-2  photostim.ulator , located 
12  in.  in  front  of  the  nasion,  was  used  to  deliver  flashes 
(10  msec  duration  at  intensity  4)  through  the  subject's 
closed  eyelids.  The  flash  tube  was  housed  in  a sound- 
attenuating chamber,  and  wiiite  noise  was  used  to  mask  any 
flash-tube  sounds  that  might  have  escaped  the  attenuator. 

In  addition,  the  subject  was  kept  both  alert  and  distracted 
from  the  flash  stimuli  by  being  engaged  in  a two-way  conver- 
sation with  a laboratory  assistant  throughout  the  experi- 
mental session. 

Digital  data  processing  in  this  study  vses  carried  out 
on  an  SEL-840  computer.  Prior  to  the  cxiseriment,  the  com- 
puter was  programmed  to  generate  a sequercr  of  inteivals 
with  values  between  2 sec  and  8 sec  selected  from  a table 
of  random  permutations.  This  variation  in  interval  size 
was  included  to  minimize  tht'  influence  of  expectancy  wb.ich 
might  develop  from  fixed  intervals  of  stimulus  {rcsentat icn 
and  to  counteract  the  possibility  that  some  lingering  stimulus- 
locked  components  in  the  ERG  m.ight  contaminate  the  phase- 
trigger  time  relationship  in  the  st  in.ulus-OFF  control  trials. 
These  time  intervals  were  recorded  on  analog  tajc  alone  with 
flash  commands  (st imulus-ON)  and  trigger  commands  ( st i nu 1 us-OFF ) 


in  separate  tape  channels.  The  beginning  of  each  interval 
had  associated  with  it  either  a stimulus  command  or  a trig- 
ger command.  Flash  and  no-flash  trials  (i.e.,  intervals) 
occurred  in  simple  alternation.  This  tapt'  was  then  used  to 
control  the  stimulus  generator  during  the  experiments  and 
to  provide  signals  representing  experimental  parameters  to 
be  recorded  onto  a second  analog  recorder  along  with  the  EEG 
signals.  These  analog  recordings  were  then  digitized  at 
2-msec  intervals  and  stored  on  digital  taj-e  to  await  further 
analysis . 

Obviously,  if  one  wishes  to  study  alpha  phase,  it  is 
of  great  importance  that  care  be  exercised  to  avoid  phase 
distortion  in  the  processing  of  the  EEC  signals.  Therefore, 
the?  digitized  EEG  data  were  processed  using  a nonrecursioe 
phase-distortionless  digital  filter.  The  parameters  of  this 
transversal  filter  were  (a)  bandpass  of  6.0  to  18.0  Hz; 

(b)  length  of  151  samples  at  2 msec  per  sample;  and  (c) 

Hamming  weighting.  The  delay  added  to  the  EEG  data  from,  the 
action  of  this  noncausal  filter  was  also  introduced  into  the 
stimulus  channel  and  the  trigger  channel  to  preserve  the 
original  time  relationship  between  stimulus  ccrmiands  (or 
trigger  commands)  and  the  EEG  tim.t?  functions. 

In  subsequent  discussions  wc?  shall  refer  to  these 
digital  data  in  terms  of  sample-epochs  or  samp  le-funct  ior^s 
corresponding  to  the  various  trials.  The  phrase  "st  im;ulus-ON 
sample-function"  will  refer  to  a sample-epoch  that  begins  with 
a stimulus  (flash)  and  is  sampled  regularly  at  2-msec  in'.erwil 


The  phrase  "stimulus-OFF  sample-function"  will  refer  to  a 
similar  epoch  that  begins  with  a trigger  only  (indicating 
that  the  stimulus  was  "OF,"  i.e.,  omitted). 

In  th  'resent  work  we  utilize  the  quadrature  defini- 
tion of  phase  as  originally  introduced  by  Weaver  (1956)  and 
applied  by  Ein-Gal  and  Lai  (1973)  to  EEG  recordings  sup- 
plied by  Anliker.  Briefly,  the  quadrature  procedure  is  as 
follows;  a narrowband  time  function  can  be  described  in 
terms  of  frequency,  amplitude,  and  phase  values.  If  the 
frequency  can  be  specified  (either  as  a constant  or  a 
tracked  variable) , it  is  possible  through  quadrature  analy- 
sis to  describe  the  instantaneous  phase  and  amplitude  of 
the  input  signal  with  reference  to  coherent  sine  and  cosine 
reference  waves  of  the  same  mean  frequency  as  the  signal. 

In  our  experiments  the  mean  alpha  frequency  was  obtained 
from  autocorrelation  analysis  of  the  EEG  data,  the  dominant 
period  of  the  correlogram  being  taken  as  the  best  estimate 
of  the  mean  alpha  period.  Quadrature  analysis,  using  this 
frequency,  compares  the  variable  signal  with  sine  and  cosine 
coherent  references  of  the  defined  frequency  and  yields  an 
in-phase  component  and  a quadrature  component.  These  two 
components  are  used  to  estimate  the  amplitude  and  the  phase 
of  the  input  sample-functions.  In  the  present  application, 
we  have  made  the  further  specification  that  the  positive 
peak  of  the  cosine  reference  wave  always  coincides  with  the 
beginning  (time  zero)  of  each  sample-function,  i.e.,  coincides 
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with  the  moment  at  which  either  a stimulus  or  a trigger- 
only  control  is  activated.  The  contingent  phase  of  a par- 
ticular sample-function  is  defined  in  terms  of  the  phase 
contingency  range,  of  which  there  are  20  equal  18°  ranges, 
within  which  the  quadrature  value  of  the  sample-function 
falls  at  time  zero.  Based  upon  this  quadrature  method  for 
the  classification  of  alpha  phase,  we  have  developed  and 
applied  two  forms  of  phase  contingent  signal  analysis  which 
we  shall  call  (1)  phase  contingent  averaging  and  (2)  phase 
contingent  expansion. 

RESULTS  AND  DISCUSSION 

It  is  instructive  to  look  first  at  the  stimulus-OFF 
(i.e.,  trigger-only)  phase  contingent  averages  in  Figs.  FI 
and  P2.  It  will  be  observed  that  the  phase  contingent 
averages  are  all  quite  substantial  nonzero  functions.  They 
are,  furthermore,  all  quite  similar  except  for  the  phase 
differences  at  time  zero.  The  peaks  of  the  anplitude  en- 
velopes of  these  phase  contingent  averages  are  all  located 
at  time  zero.  This  is  to  be  expected  because  these  phase 
contingent  averages  are  equivalent  to  the  triggered  averages 
of  sato  (1962)  and  the  triggered  correlations  of  Boer 

and  Kuyper  (1968).  Consequently,  these  phase  cotitingcnt 
averages  partake  of  certain  basic  foature-s  of  the  auto- 
correlogram.  For  instance,  there  is  a characteristic  fall- 
off  in  the  amplitude  with  increasing  time  distance  from, 
time  zero.  This  is  accounted  for  in  our  phase-incoherence 
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SUBJECT:  K.L. 
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Phase-  continnent  expansions  ol  AVEP  data  i i-oni  subieet  K.I..  The  expansion  -I 
St  i mu  Uis-OFF  (ti-it'i;ei'  but  no  stimulus)  data  is  on  the  loTt  side  ol  the  1 i t^u  re  , 
the  st  iinulus-ON  expansion,  on  the  rii;hl  side.  The  various  contint;eiu  .ive  r:ii;es 
are  stacked  as  a vertical  sei‘ies,  with  continKencv  I on  the  bottom  and  eon- 
tinBoncy  2(J  at  the  top.  The  Grand  Mean,  or  mean  of  the  eonl  indent  t ans, 
equivalent  t.t  the  noncon  t i iiBent  .WTIP,  is  located  immediately  aliove  i 'in!  i I4;oik  \ 
lin  and  Is  iiuhealed  liy  a heaviei-  line.  The  stimulus  ('M-  iriBB'-r)  tire  is  -.lioun 
.It  the  top.  iTie  averaKi:  numbi-r  ot  sam|)U-- lunet  ions  in  each  eont  indent  .m  iai-.e 
IS  1(1. 
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HGURK  I-.2. 

Phase  continKOiit  expansions  oi  AVEP  data  i rom  subjec  t 1..U.  The  expansion  ol 
stiniulus-OFF  (trigger  but  no  stimulus)  data  is  on  the  li'ft  side  o(  the  ligure 
the  stimulus-ON  expansion,  on  the  right  side.  The  various  coiUingen;  avc'rages 
are  stacked  as  a vertical  series  beginning  with  contingency  1 on  the  bot t . 

The  Grand  Mean,  or  mean  ol  the  contingent  means,  equivalent  to  the  noncontiti 
gent  AVEP,  is  located  immediately  above  contingency  20  and  is  indicated  Ijy  a 
heavier  line.  The  stimulus  (or  ti’igger)  time  is  shown  al  t lie  top.  I'he  .ivi-rage 
number  ol  sample-liinct  ions  in  each  contingent  avc'rage  is  to. 


model  (Floyd  ^ 1973)  by  the  propagation  of  phase  in- 

coherence as  a function  of  increasing  time  distance  from 
the  classification  point  which  is  located  at  time  zero. 

For  a coherent  wave  train  there  is  no  fall-off  in  the  ampli- 
tude of  the  envelope  of  phase-triggered  tim.e  averages  (or 
of  the  envelope  of  the  autocorrelation  function) . Alpha 
may  be  described  as  a quasi-coherent  rhythm.,  which  means, 
according  to  the  phase-incoherence  model,  that  the  rhythm; 
has  a coherent  mean  frequency  and  a limited  amount  of  phase 
variance  or  phase  incoherence;  the  greater  the  observed 
phase  variance,  the  faster  the  rate  of  fall-off  in  the  ampli- 
tude envelope  of  the  phase  contingent  average.  It  will  be 
observed  in  these  phase  contingent,  stimulus-OFF  averages 
that  there  is  sufficient  coherence  to  "predict"  some  rather 
definite  phase  values  at  those  latencies  in  which  the  m.ajor 
evoked  deflections  appear  in  the  AVER  (which  is  equivalent 
to  the  Grand  Mean  of  the  phase  contingent  averages)  . V.'hon 
these  phase  "predictions"  in  the  stimulus-OFF  contingent 
averages  are  compared  with  the  phase  properties  exhibited 
by  their  corresponding  stimulus-ON  contingent  averages, 
there  are  some  striking  changes.  To  facilitate  visual  exami- 
nation  of  these  rather  complex  contingent  expansions.,  we 
have  provided,  in  Figs.  P3  and  P4,  the  same  data  as  in  Figs. 
PI  and  P2;  some  schematic  lines  have  been  added  to  call 
attention  to  those  features  that  appear  to  us  to  be  moved 


as  a result  of  stimulation.  We  have  used  dotted  lines  to 
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iigurf:  I’.si. 

ScliemnLic  ovorlays  on  tliu  data  ol  Hubjuct  K.I..  ( i ii  I'lt;.  I’.l).  riu'sc  ' > vi- r 1 a \ :. 

indicate  t hC!  principal  d i 1 1 c iimccs  to  Ik;  oliscrvcd  between  t tu‘  a I i I'lu  1 n s-Ol'I- 
and  t ho  stimuJus-ON  conditions.  Notice  the  phase  sli  i I I ol  the  Nl-Pii-Nti  le 
lioncnts  If  om  a diagonal  oric-nlation  in  t lie  s I i mu  1 us-on-  expansion  to  a ve  r' 
it:al  oricMilation  in  l In;  stiniuUis-ON  expansion,  .oakiiip  'he.'-.e  loi'ponent  . in 
all  the  coni  1 nnoni' 1 os  ap;ree  mori'  or  less  with  the  Grand  Mean. 
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SUBJECT:  L.H, 


I I'lciuiU':  p.i. 

[ Schematic  overlays  on  Ihe  coni  ^n^;l;nt  ex|>;insions  ol  subject  b. II.  (lit;.  I’.J). 

[ These  overltiys  Imlictile  the  principal  d i I le  renc-cs  to  be  i ibsc' rt  c-d  bclwi’cn 

j the  St  i mu  1 us-OKF  and  the  stimulus-ON  conditions.  Notice  the  jihase  sh  i 1 t .■! 

I the  N1-P2  N2  cotnixmenls  i rom  a diagonal  ■>  I'l  cut  a t i on  in  the  st  i imi  I us-oi'T 

expansion  (leitl  to  a vertical  orientation  in  the  st  i mi  liis-ON  expansion 
(right),  making  these  comiionents  in  all  contingencies  agree  n.o  re  r h . 
with  Ihe  Grand  Mi'an.  Holh  subjects  show  essentially  t lie  ,aiie  e||e<i  ■ | 
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indicate  the  crests  or  surface-negative  peaks  and  solid  lines 
to  indie rf  the  troughs  or  surface-positive  peaks.  In  the 
case  of  the  stimulus-OFF  contingent  averages,  it  will  be 
notices  that  these  schematic  lines  are  all  essentially  paral- 
lel diagonals  extending  uj^ward  and  to  the  left  from  contin- 
gency 1 to  contingency  20  and  connecting  corresponding  phase 
points  en  route.  However,  when  we  examine  the  stimulus-ON 
phase  contingent  averages,  we  notice  that  these  lines  are 
moved  in  such  a manner  that  the  Nl,  P2,  and  N2  lines  are 
swung  into  more  or  less  vertical  orientation  under  the  Nl, 

P2,  and  N2  deflections  in  the  Grand  Mean  (i.e.,  AVEP)  of 
the  contingent  averages.  You  will  notice  that  the  effects 
are  quite  similar  in  both  subjects.  This  surprisingly 
systematic  alteration  of  the  phase  properties  is  inter- 
preted as  strong  support  for  the  view  that  the  photic  stimuli 
act  to  produce  a phase  modulation  of  the  alpha  rhythm  toward 
a preferred  phase  which  is  consistent  with  the  major  Ni-P2-N'2 
deflections  in  the  AVEP.  In  other  words,  it  is  this  phase 
shifting  of  the  alpha  from  the  autonomous-phase  locations 
to  the  evoked-phase  locations  which  accounts  for  this  por- 
tion of  the  AVEP.  That  is,  the  "signal"  or  evoked  component 
is  nothing  more  than  the  phase  sliifting  of  the  autonomous 
phases  to  the  evoked  phases.  If  this  interpretation  is 
correct,  there  is  no  need  to  look  for  a nonalpha  "sicjnal" 
to  account  for  this  porti<5n  of  the  AVKl' . In  an  earlier 
publication,  Floyd  et  a_l_.  (1^73)  presented  data  in  support 

of  the  idea  that  the  so-called  " a f t or-d  i srha  rtie"  or 
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"reverberation"  which  is  frequently  observed  at  longer 
latencies  could  be  accounted  for  solely  on  the  basis  of 
phase  locking  which  had  occurred  earlier  in  the  evoked 
response,  i.e.,  at  the  time  of  the  appearance  of  the  N1-P2-M2 
deflections.  In  the  present  work  we  have  produced  evidence 
that  such  phase  alignment  does  actually  occur  in  this  region 
of  the  AVEP. 

From  their  study  of  auditory  evoked  potentials  using 
amplitude-spectral  and  phase-spectral  methods,  Sayers  e;t  al. 
(1974)  concluded  that  "a  moderate-level  stimulus  that  is 
effective  in  evoking  an  evident  response  apparently  does  not 
do  so  by  contributing  an  obvious  additive  component."  And 
they  go  on  to  say  that  "an  effective  stim.ulus  acts  m.ainly 
to  phase-control  the  spontaneous  activity."  Their  view  is 
somewhat  similar  to  ours  although  they  do  not  identify  the 
spontaneous  activity  as  alpha.  Also,  whereas  our  model  is 
based  upon  a single  frequency  component  (alpi)ia)  , their  model 
seems  to  require  the  existence  of  a whole  series  of  coherent 
oscillators  (at  least  10  harmonics),  each  with  the  capacity 
for  being  individually  adjusted  in  phase  and  amplitude.  It 
will  be  most  interesting  to  see  wfu>ther  the  auditory  evoked 
response  can  be  successfully  analyzed  in  termts  of  phase  modu- 
lation of  a single  frequenev'  component  su<-h  .i.s  alpha.  Tlie 
findings  of  Sayers  et  c^.  (1974)  certainly  suggest  that  the 

matter  is  worthy  of  careful  inv'est  i gat  ion . It  is  tempting 
to  speculate  that  analogous  phase  encoding  mechanisms  ray 
operate  in  all  of  the  cerebral  sensory  me'dal  i t ies . 
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AN  ALPHA  PHASE  MODULATION/DEMODULATION  THEORY  OF  CEREBRAL 
ENC:ODING/DEC 

In  contrast  to  the  widespread  assumption  that  alpha 
is  nothing  but  noise,  we  (Anliker  and  Floyd,  1977a)  have 
proposed  a radically  different  concept  whicli  confers  upon 
alpha  activities  a central  role  in  the  encoding  and  de- 
coding of  sensory  information.  We  call  our  concept  the 
phase  modulation  theory  of  cerebral  encoding  and  decoding. 
Accordi  g to  this  theory,  alpha  serves  as  a phase  coherent 
carrier  on  which  the  stimulus  effect  is  encoded  as  a phase 
modulation.  That  is,  the  autonomous  alpha  represents  the 
carrier  without  modulation  while  the  stimulus-influenced 
alpha  represents  the  carrier  with  sensory  modulation.  The 
autonomous  cortical  alpha  is  considered  to  represent  the 
cortical-following  response  to  the  thalamic  pacemaker  (the 
latter  being  buffe''ed  fron  specific  sensory  influences). 
However,  the  cortical  following  of  the  thalamic  pacem.aker 
is  modulat<‘d  in  the  presence  of  appropriate  sensory  input  to 
the  cortex.  The  phase  modulatioi\  theory  states  that  through 
som'  process  of  comparing  the  pacemaker  signal  to  the  res- 
ponse of  the  cor  respond ina  sensory  cortex,  the  monitors  (pre- 
sum.ably  corticcjl)  are  able  to  distinciuish  between  (a)  endo- 
genous activities,  i.e.,  those  cortical  responses  which  can. 
be  .iccounted  for  solely  on  the  basis  of  pacemaker  following, 
and  (b)  exogenous  activities,  i.e.,  those  cortical  responses 
which  differ  from  the  pacemakt  i i nst  ructi  oi.s . According  to 
the  phase  modulation  theory,  the  phase-ivodul  ated  component  s 


are  culled  out  for  further  processing  while  the  endogenous 
or  unmodulated  phase  components  are  disregarded.  In  other 
words,  the  endogenous  activity  serves  as  the  undifferentiated 
medium  upon  which  external  influences  are  encoded  as  phase 
shifts  or  modulations  of  the  autonomous  phase.  This  is  re- 
garded as  a crucial  step  in  the  incorporation  of  exogenous 
influences  into  cerebral  processes.  One  might  say  that  the 
autonomous  activities  constitute  a low  priority  input  to 
the  cerebral  incorporation  process  and  are  ordinarily  dis- 
regarded (i.e.,  nulled  or  filtered  out)  whereas  the  exo- 
genous activities  constitute  high  priority  inputs  which  will 
tend  to  be  accepted  (passt^d  upward  in  the  neural  hierarchy) 
for  higher  level  processing  in  attention  and  memory.  How 
many  filters  the  input  must  pass  before  emerging  before  the 
"throne  of  consciousness"  is  a question  which  the  phase 
modulation  theory  does  not  attempt  to  answer.  But  the  theory 
does  state  that  the  phase  modulation  is  the  sine  qua  non  of 
cerebral  incorporation,  tlie  critical  step  wherein  the  exo- 
genous influences  are  translated  into  a code  or  format 
which  is  the  basic  language  of  the  cerebrum.  This  is  not  to 
say  that  frequency  and  am)-ilitude  influences  ar(>  ignored  ly 
the  cerebrum  (for  obviously  they  are'  not):  however,  it  does 
mean  that  the  cerebral  "carrier"  is  pacemaker  phas(>  and  all 
time-domain  inputs  are  encoded  as  phase  m.odul  a t i ci'. . An.othei- 
prediction  of  the  phase  modulation  theory  is  that  the  de- 
coding of  the  stored  phase-modulated  signals  will  reejuire  the 
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playback  will  require  the  reactivation  of  the  carrier  (i.e., 
the  presence  of  the  autonomous  alpha  rhythm)  for  the  re- 
awakening of  mnemonic  records,  i.e.,  for  remembrance.  The 
phase  modulation  theory  predicts  that  rememi^rance  will  be' 
distorted  unless  the  carrier  frequency  during  playback  is 
identical  to  that  during  recordina.  This  would  account  for 
many  of  the  perceptual  distortions  which  are  observed  when 
the  autonomous  alpha  frequency  deviates  significantly  fresm 
its  normal  alert  value.  Such  deviations  are  observed  in 
normal  individuals  during  states  of  drowsiness  and  dream- 
ing, and  they  can  result  from  the  influence  of  drugs,  trau- 
matic brain  damage,  and  metabolic  disorders.  It  should  be 
noted  that  it  is  the  normal  alert  alpha  that  exhibits  the 
greatest  phase  coherence;  when  alpha  departs  from  its  normal 
alert  frequency  there  is  a reduction  in  the  phase  coherence 
of  the  autonomous  rhythm.  Although  there  are  no  conclusive 
EEG  records  for  the  alpha  frequency  differences  in  psychiatric 
patients  during  normal  periods  of  contact  and  during  periods 
of  disorientation,  the  theory  would  predict  a fre<;uency 
difference.  The  low  alpha  index  that  characterizes  schizo- 
phrenics may  be  a clue  that  the  norm.al  relationship  between 
the  pacemaker  and  the  cortex  is  disturbed;  if  so,  this  would 
tend  to  cause  distortions  in  the  perceptual  processes. 

It  seems  to  us,  however,  that  some  cautionary  notes 
are  in  order.  The  alpha  phase  modulation  theory  in  its 
present  form  does  not  postulate  a single  alpha  paccm.akei 
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for  the  entire  cortex.  Rather,  we  postulate  that  each 
portion  of  the  cortex  has  its  own  thalamic  pacemaker  refer- 
ence along  lines  suggested  by  Andersen  and  Andersson  (1968). 
We  would  also  point  out  that  our  model  in  its  present  form 
is  concerned  primarily  with  the  time  domain;  we  have  not  at 
this  time  attempted  to  reconcile  our  theory  with  an  essen- 
tially spatial  model  such  as  the  holographic  model  of  cere- 
bral encoding  proposed  by  Pribram  (1971)  except  to  observe 
that  the  quasi-coherent  alpha  activity  could  provide  the 
coherent  reference  which  is  required  by  a holographic  systcmi. 
Finally,  we  believe  that  a careful  distinction  must  be  made 
between  the  use  of  the  term  "phase"  in  referring  to  the 
encoding  of  the  sensory  data  of  visual  space  (see  Pollen 
et  al . 1971) , and  our  use  of  "phase"  as  a parameter  of  the 
time  domain.  In  contrast,  having  already  utilized  the  time 
domain  to  account  for  the  encoding  of  spatial  phase  infor- 
mation (Pollen  e^  1971),  Pollen  and  Trachtenberg  (1972) 

are  attracted  to  the  alpha-as-noise  hypothesis.  They  state 
that  "even  if  it  were  assumed  that  all  colls  in  the  visual 
cortex  were  rhythmically  excited  by  alpha  activity  and  in 
phase  with  each  other,  simple  cell  responses  to  visual 
stimuli  in  different  parts  of  the  receptive  field  would 
arrive  at  the  complex  cells  with  different  latencies  and 
thus  at  different  periods  of  the  alpha  cycle.  The  alpha 
activity  would  facilitate  some  inputs  and  inhibit  others 
in  a nonpr(?dictable  way--a  noise  process.  Thus  alpha  block, 
wht'ther  by  a central  active  inhibitory  mecl.anism  c'r  by  a 


more  general  desynchronization  of  rhythmic  activity,  might 
serve  to  reduce  a neural  noise  level."  We  believe  that  our 
phase  modulation  theory  provides  a viable  alternative  ex- 
planation to  that  advanced  by  Pollen  and  Trachtenberg  (1972) 
for  the  alpha  blockage  observed  during  visual  attention  to 
nonuniform  surfaces.  That  is,  the  desynchronization  of  the 
cortical  response  may  represent  the  complex  phase  modula- 
tion of  a large  cortical  area  in  response  to  complex  spatial 
patterns  of  sensory  input.  These  diverse  phase  adjustments 
of  the  cellular  "alpha"  activities  could  be  expected  to 
greatly  increase  the  phase  variance  of  the  cell  pool  and 
thereby  reduce  the  m.agnitude  of  the  potentials  measured  at 
a scalp  electrode.  By  contrast,  in  the  absence  of  sensory- 
phase  modulation,  a much  greater  phase  coherence  in  the 
cortical  activity  is  expected  because  of  coherence  of  the 
thalamic  pacemaker. 

We  have  concluded,  therefore,  that  sirriple  timie  acjeraqing 
of  sample  functions  of  EF.G  data  confounds  the  influence  of 
alpha  phase  contingency  in  the  generation  of  the  visual  evoked 
response.  By  classifying  each  sample-function  according  to 
the  phase  of  the  alpha  rhythm  at  the  moment  of  stimulus  de- 
livery, it  is  possible  to  regroup  the  set  of  sample-f unct ions 
into  subsets  representing  various  phase  contingency  ranges. 

The  time  averages  of  these  subsets  constitute  a phase  con- 
tingent expansion  of  AVEP . The  results  of  this  analysis  can 
be  interpreted  as  showinc:  that  tlu'  visual  stimulus  (flash) 
acts  to  phase-mcxlulate  tlie  autonomous  alpha  rhythm  and  that 


this  effect  accounts  for  the  principal  N1-P2-N2  component 
of  the  AVEP. 


ALPHA  PHASE  PROBABILITY  ANALYSIS  OF  THE  AVERAGE  VISUAL 
EVOKED  RESPONSE 

The  present  study  was  designed  to  extend  the  alpha 
phase  contingency  analysis  and  to  produce  a more  quantita- 
tiv€i  evaluation  of  the  phase  contingency  effects  by  re- 
analyzing AVER  data  from  the  standpoint  of  phase  contingent 
probabilities.  We  believe  that  the  alpha  phase  contingency 
classificatory  methods  described  here  represent  a definite 
advc'ince  in  the  techniques  available  for  the  analysis  of 
evoked  responses  and  for  the  quantitative  characterization 
of  spontaneous  alpha  waves. 

Phase -Pis tort ion less  Filtering 

In  the  analysis  of  alpha  phase  it  is  a matter  of  ut- 
most importance  to  take  suitable  precautions  to  avoid  dis- 
torting phase  in  the  process  of  filterino  the  Er,G  signals. 


Therefore,  in  the  present  study  the  digitized  EEG  data  were 
processed  using  a nonrecursive  phase-distortionless  filter. 
The  parameters  of  this  transversal  filter  were  (a)  bandpass 
of  6.0  to  18.0  Hz;  (b)  length  of  151  samples  at  2 msec  per 
sample;  and  (c)  Hamming  weighting.  The  delay  added  to  the 
EEC  data  from  the  action  of  the  noncausal  filter  was  also 
introduced  into  the  stimulus  record  and  into  the  trigger 
record  so  as  to  preserve  the  original  time  relationship  be- 
tween stimulus-commands  (or  trigger-commands)  and  the  EEG 
sample  functions. 

DATA-PROCESSING  METHODS 

Quadrature  Analysis 

In  this  study  we  have  employed  the  quadrature  defini- 
tion of  phase  which  can  be  described  briefly  as  follows:  a 
narrowband  time  function  can  be  analyzed  into  frequency,  ampl 
tude,  and  phase  variables.  If  the  frequency  can  be  speci- 
fied (either  as  a constant  or  as  a trackable  variable),  it 
is  possible,  through  quadrature  analysis,  to  describe  tht 
instantaneous  phase  and  amplitude  of  the  input  signal  with 
reference  to  coherent  sine  and  cosine  waves  of  the  same  mean 
frequency  as  the  signal.  In  our  study,  the  mean  alpha  fre- 
quency was  obtained  by  autocorrelation  analysis  of  the 
stimulus-OFF  sample  functions  of  EFG , the  dominant  period 
of  this  autocorrelogram  being  taken  as  the  best  estimate  of 
the  mean  alpha  period.  Quadrature  analysis,  using  this  fre- 
quency, compares  the  variable  input  si  anal  with  sine-  and 


cosine  references  of  the  estimated  frequency  and  yields  an 
in-phase  and  a quadrature  component.  These  two  components 
are  used  to  estimate  the  amplitude  and  the  phase  of  the 
input  sample  functions.  In  the  present  application,  we  have 
made  the  further  specification  that  the  positive  peak  of 
the  cosine  reference  wave  always  concides  with  the  begin- 
ning (time  zero)  of  each  sample  function,  i.e.,  the  time 
at  which  either  a stimulus  or  a trigger  is  activated.  The 
contingent  phase  of  each  sample  function  is  defined  in  terms 
of  the  phase  contingency  range  (we  recognized  twenty  equal 
18°  phase  ranges)  within  which  the  quadrature  phase  of  the 
sample  function  falls  at  time  zero.  A quadrature  phase  time 
function  describes  the  phase  difference  between  the  phase  of 
the  coherent  reference  wave  and  the  phase  of  the  EEG  sample 
function  at  each  post-stimulus  (or  post-trigger)  latency. 

A quadrature-phase  time  function  is  obtained  for  each 
EEG  sample  function.  This  phase  function  is  used  as  the 
basis  for  further  analyses  of  alpha  phase  properties. 

Phase  Probability  Analysis 

The  reader,  no  doubt,  is  generally  familiar  with  the 
procedures  for  constructing  frequency  (of  events)  histogram.s 
and  for  converting  these  histograms  into  phase  probability 
distributions.  The  GRA.ND  PHA.SE  PROBABILITY  MATRIX  contained 
in  each  of  the  four  figures  (part  E)  is  really  nothina  m.ore 
than  a series  of  phase  probability  distributions  in  which 
each  of  the  probability  distributions  desciibes  a set  of 
quadrature  phase  values  observed  at  a particular  latency.  In 


the  absence  of  stimulation  one  would  expect  that  a suffi- 
ciently large  population  of  quadrature  phase  functions  would 
generate  a phase  probability  matrix  with  an  approximately 
flat  surface  because  all  phases  of  alpha  should  be  equally 
probable  at  every  latency.  Any  significant  ard  consistent 
deformation  of  this  surface  would  be  considered  as  evidence 
for  the  existence  of  some  sort  of  phase  constraint. 

Inasmuch  as  the  grand  phase  probability  matrix  provides 
no  account  of  phase  contingency  effects,  it  may  preve  useful 
to  construct  phase  contingent  phase  probabi] ity  matrices 
as  well.  Since  our  phase  classification  scheme  recognizes 
20  different  phase  contingency  ranges,  it  is  feasible  to 
produce  a separate  phase  probability  matrix  for  each  of 
these  contingencies.  v;e  have,  in  fact,  carried  out  this 
analysis.  What  is  the  general  shape  of  the  phase  contingent 
phase  probability  matrix?  In  each  case  we  should  expect  to 
obtain  a delta  function  (Dirac)  at  time  zero  because  all  of 
the  samples  would  be  within  the  same  phase  range  at  the  time 
of  classification,  by  definition.  With  increasing  tim.e 
distance  (latency)  from  the  classif ication  point,  the  propa- 
gation of  phase  incoherence  would  cause  the  probability  dis- 
tributions at  longer  latencies  to  be  less  peal^ed  than  at 
earlier  latencies  because  phase  variance  increases  as  a 
function  of  latency  magnitude.  On  the  other  hand,  in  the 
absence  of  stimulus-induced  phase  effects,  this  " feigning 
out"  of  the  phase  contingent  phase  probability  matrix  should 
be  orderly  and  symmetric^ll . Consequently,  any  consistent 
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and  significant  differences  observed  between  the  correspond- 
ing stimulus-ON  and  stimulus-OFF  phase  contingent  phase 
probability  matrices  could  be  interpreted  as  evidence  for 
the  existence  of  a stimulus-induced  phase  probability  change. 

The  grand  phase  probability  matrices  in  Figs.Cl-4  are 
displayed  in  pseudo  three-dimensional  form.  While  this  is 
attractive  to  look  at,  this  display  distorts  probability 
along  the  latency  axis.  To  overcome  this  limitation,  two 
additional  plots  are  provided:  modal  phase  (part  C)  which 
shows  the  location  of  the  probability  peaks  when  projected 
onto  the  phase-latency  surface;  peak  amplitude  (part  B) 
which  shows  the  amplitude  (relative)  of  the  probability 
peaks  above  the  average  probability  surface. 

It  wi'l  be  helpful  to  review  briefly  the  m.anner  in 
which  the  COMPOSITE  PLOT  (part  D)  of  the  model  (cr  peak) 
phases  of  twenty  (a  complete  set)  phase  contingent  phase 
probability  matrices.  Since  it  is  not  practical  to  display 
80  phase  contingent  phase  probability  plots,  these  data  are' 
summarized  by  extracting  the  model  phase  function  (i.e.,  the 
peak  phase  at  each  latency)  from  each  of  the  matrices  and 
entering  these  points  into  a common  phase-) atency  display. 

By  comparing  the  composite  plot  with  its  associated  orand 
phase  probability  matrix  it  is  possible  to  discern  how  the 
various  contingent  functions  respond  to  the  stimulus.  These 
effects  will  be  mentioned  in  the  Results  section. 
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RESULTS 

Our  analyses  of  the?  experimental  data  are  displayed 
in  Figs.  Cl-4.  These  figures  are  designed  to  help  the  reader 
make  his  own  comparisons  of  the  various  analyses  contained 
within  each  figure  and  also  to  aid  his  comparisons  of  the 
corresponding  portions  of  the  different  figures.  Obviously, 
there  are  many  interesting  aspects  of  these  data.  Yet,  it 
will  be  necessary  to  confine  our  attention  to  the  most 
salient  features  of  these  data. 

Figures  Cl  and  C2  contain  the  stimulus-OFF  and  stimulus- 
ON  data  of  subject  K,  L.  The  AVER  in  Fig.  Cla  is  close  to 
zero  amplitude  throughout,  which  is  to  be  expected  of  stimulus- 
OFF  data.  Notice  that  all  of  the  other  plots  are  consistent 
with  this.  The  peak  amplitude  plot  (B)  shows  that  the 
largest  non-zero  value  in  the  grand  phase  probability  matrix 
is  quite  small.  The  composite  plot  (D)  shows  that  the  con- 
tingent modes  are  not  sicini  f icantly  phase  constrained.  how- 
ever, we  see  quite  a different  pictuie  when  we  exc»mine  Fig. 

C2  (stimulus-ON)  . There  is  a substantial  Nl-}''2-N2  response 
in  the  AVER  (A).  Tn  the  grand  pliase  probability  jjlot  (E) 
there  is  a large  deformat.ion  of  the  surface  of  the  matrix  in 
the  latency  domain  that  corresponds  with  the  evoked  response. 
The  magnitude  of  this  response  can  be  appreciated  by  looking 
at  the  peak  amplitude  plot  (B)  . Wliat  is  of  great  intcrtst 
to  us  is  the  convergence  of  ttie  modal  phast^s  of  the  con- 
'ingent  probability  matrices  (D)  at  late'ncics  correspond ino 
'■  the  evokf'd  response  (!Jl-P2-N2).  This  stromjly  suggests 

eg 


TRIGGfcR 


STIMULUS  OFF 


TIME  IN  4U  INTERVALS 


AVERAGE  VISUAL  EVOKED  RESPONSE 


PEAK  AMPLITUDE  OE  GRAND  PHASE  PROBABILITY  MATRIX 


MODAL  PHASE  OF  GRAND  PHASE  PROBABILIIY  MATRIX 


TIME  4 mvK  INTERVALS 


COMPOSITE  PLOT  MODAL  PHASES  OF  TWENTY  CONTINGENT  PHASE  PROBABILITY  MATRICES 


GRAND  PHASt  PROBABILITY  MATRIX 


^ Tv(„v  , ; ! •F'  '''Av  ■ 

180  ■_.>  '• ; 

0 40  80  l^'O  Tw  Too  240  2H0  .^20  ~360~  400  440  480  h20  ‘.RO  WK’  840  880  ’20 

»>>MN 


I’KHIRK  C' . 1 . A I'll  I 1 ixa  I I lit!  nl  I hi'  vai'ums  aii.ilysi'S  ■'!  •iiliiiil  K.I.'- 
st  1 mil  1 us-OKF  (kila  nvi'P  a cuTiiiniin  t i iiu'-hasi' . K.iiTi  nl  tln'  analyst's  is  nli'iit 

lii'd  hy  a capital  li'llcr  (A  - Kl  liKali'd  in  thi'  lilt  ' si  ri;  i n .iiul  li\  a 

ill 'SC  r 1 pt  i vi'  label  placi'il  i mmi'il  i a t c I v .ihnvc  Ihc  i;r.ii)liic  i.ateiial  . F'  i 

-'xaiiipli',  A AVKIi-AtiK  VI SUAI  KVOKKD  KKSlx )NS!' 


STIMULUS  ON 


TIME  IN  40  INUHVALS 


AVERAGE  VISUAL  EVOKCU  RESPONSE 


A 0^  - 


PEAK  AMPLITUDE  OF  PROBABILITY  MATRIX 


MODAL  PHASE  OF  PROBABILITY  MATRIX 


r 5 0. 

I 360 


riME  4 mvet  INTEIWALS 


COMPOSITE  PLOT  MODAL  PHASES  OF  I WENTY  CONTINCU  NT  PHASE  PROBABILITV  MATRICES 


10 


n ^ ’ ~ 


* ’V 


GRAND  PHASE  PROt^AHILITV  MATRIX 

180  • rr-F 

i ■ Jf  ' iT  , jf^ 


E I TWi 


//•/  . ' 





I)  41)  HO  l?n  160  /OH  /4I1  /Hll  I/O  IMl  4()(l  44|i  4BM  'p/0  ‘■60  PplHI  1.4"  6HI'  .'■• 

IIMI  ■■■>■■■ 


lIGUlU','  C.2.  CO  1 1 (k;i  ! I I'll  "I  ttio  v.irioiis  jiikiIvsi";  'I  .iiih)i'<-l  k.I,.'- 
slimulus-ON  ilatii  over  a ii  imiin  ui  I i iiio  - lia  ■;<-• . k.icli  "I  tho  .iiialvsi's  i,  i (Iimi  1 i I 1 1 il 
liy  a capital  Irtlci  (A  - K>  locati'd  in  I lio  Iclt  ".iii'.in  ami  hy  .i  ili-  a i i pt  i m- 
label  place'll  i mini'd  i a I '■  1 y alvvc  the  ^^rapliic  na  1 1 ' i- 1 a I . 1 i i\.i':ipli' 

K (iltANII  PIIASK  PKOllAIll  I.ITY  \WTHI\. 


71 


STIMULUS  OFF 


trigger 


SUKJfcCT  t M 


TIME  IN  40  m*«c  INTERVALS 


average  visual  evoked  response 


PEAK  AMPLITUDE  OF  GRAND  PHASE  PROBABILITY  MATRIX 


MODAL  PHASE  OF  PROBABILITY  MATRIX 


I IML  a uit.  I^  I [IIVAIS 


COMPOSITE  PLOT  MODAL  PHASES  OF  TWENTY  CONTINGENT  PHASE  PROBABILITY  MATRICES 


GRAND  PHASE  PROBABILITY  MATRIX 


FIGUKK  C.;i.  A L-i)  1 1 L loiv  ul  t ho  various  analyses  n| 
sL  i mu  lus-OFK  data  over  a rnminon  t i im'-hasc  . K.uh  i>I  tin 
I tod  by  a ca()il,al  lottor  (A  - K)  liualoil  in  itu'  li’ll 
doscrl  pi  1 VO  labi'l  plarod  i niiiuMi  i a 1 o V y ah.ivo  oai  h I vpo  ' 
oxamplo.  ('  \K)I)AI.  I’tlASK  PIlOllAllTMrY  UM'KIX. 


STIMULUS  ON 


STJMULUS 


SUBJtCT  L H 


TiMt  1^  40  mwt  INTERVALS 


AVERAGE  VISUAL  EVOKED  RESPONSE 


PEAK  AMPLITUDE  OF  GRAND  PHASE  PROBABILITY  MATRIX 


MODAL  PHASE  OF  GRAND  PHASE  PROBABIl  ITY  MATRIX 


TIME  4 in$tH  INTE  HV  AL  S 


COMPOSITE  PLOT  MODAL  PHASESOF  TWENTY  CONTINGENT  PHASE  PROBABILITY  MATRICES 


GRAND  PHASl  PROBABILITY  MATRIX 


KlCiUKl-'.  C’ . 1 . A i-u  1 li)t;i  I i on  ovor  n tinmiion  I i mo-b.i  so  o)  tin-  Miiions  iiii.ilsso 
ol  subjoct  b.ll.'s  st  i mu  1 US-ON  dalu.  Kacb  o|  l lii;  aiialv'.is  is  idrnlilnd  b\ 
a oapilal  Il'ILi.t  (A  - K)  locatod  in  Uio  loll  mai'flin  and  b\  a do  a ri  pi  i M' 
label  plaood  i inniod  i a t c 1 y abov»'  each  I yjx'  ol  analysis.  Vor  ox.i'  plo,  |i 
\K)I)A1,  I’HASKS  OK  T’VVKNTV  CONTINGKNT  I’llA.SK  I’UOIIAHII.ITA-  \WTim!;;-; 


that  the  stimulus  acts  tc  shift  or  modulate  the  autonomous 


phase  toward  a preferred  phase.  j 

Examination  of  Figs.  C3  and  C4  shows  that  the  second 
subject  (L.  H.)  responded  in  very  much  the  same  manner. 

That  is,  both  subjects  show  the  convergence  of  the  modal  ; 

phases  (D)  and  the  corresponding  deformation  of  the  grand  i 

phase  probability  matrix. 

DISCUSSION 

It  seems  to  us  that  the  results  we  have  obtained  in 
these  two  subjects  are  rather  robust.  In  other  words,  there 
is  little  question  that  there  are  substantial  phase  shifts 
evident  in  the  stimulus-CN  contingent  averages;  these 
phase  adjustments  are  imposed  upon  the  avitonoir-ous  alpha 
activity,  the  latter  being  estimated  by  the  corresponding 
phase  contingent  stimulus-OFF  averages.  It  is  important  to 
take  notice  of  the  fact  that  all  of  the  phase  contingent 
stimulus-OFF  averages  are,  as  should  be  expected,  substantial 
non-zero  averages  with  a strong  alpha-frequency  rhythm.  For 
this  reason  there  is  a predicted  phase  for  the  contingent 
stimulus-OFF  (i.e.,  autonomous)  alpha  waves  iit  each  latency 
in  the  post-trigger  time  domain) . The  stimulus  appears  to 
shift  the  alpha  phase  into  conformity  with  the  evokt d Nl-r2-K2 
complex  in  the  AVER.  There  is  some  possibility  that  the 
various  alpha  phase  contingencies  contribute  differently  to 
the  various  latencies;  however,  this  possibility  will  have 
to  bo  assessed  by  further  analysis  of  the  data  in  a subse- 
quent publication.  The  reader's  attention  is  dir»ct<d  tf 
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the  earlier  remarks  conce.rning  the  effects  of  the  propaga- 
tion of  phase  incoherenc€^  on  the  shape  of  the  phase  proba- 
bility distributions  at  the  longer  latencies;  this  flatten- 
ing of  the  distribution  means  that  the  predictability  of 
alpha  phase  becomes  less  certain  at  latencies  more  distant 
from  the  classification  point.  The  iilpha  generator  is  only 
quasi-coherent . Nevertheless,  during  the  latencies  occu- 
pied by  the  N1-P2-N2  complex  of  the  AVER,  the  predictability 
of  autonomous  phase  contingent  alpha  phase  is  still  quite 
good.  That  is,  the  phase  probability  matrix  is  still  ex- 
hibiting reasonable  peakedness  at  those  latencies.  Conse- 
quently, the  bending  of  the  modal  phase  probability  func- 
tions to  produce  the  convergence  seen  in  the  stimulus-ON 
composite  plots  (part  D)  shows  the  ability  of  stimuli  to 
shift  alpha  phase  toward  a "preferred"  phase. 

In  its  present  form  our  alpha  phase  modulation  hypo- 
thesis does  not  assume  that  there  is  a single  alpha  pace- 
maker for  the  entire  cort.ex.  Instead  we  postulate  that 
each  portion  of  the  cortt;x  has  its  own  thalamic  pacemaker 
as  suggested  by  Andersen  and  Andersson  (1968).  However, 
contrary  to  their  model  v'hich  is  based  upon  barbiturate 
spindling  activity,  we  assume  that  the  waking  brain  is  charac- 
terized by  narrow  frequency  tuning,  greater  phase  coherence, 
and  more  coordinated  pacemaker  activities.  We  must  also 
point  out  that  our  alpha  phase  modulation  model  is  concerned 
primarily  with  the  time  domain.  As  yet  we  have  not  attem.pted 
to  reconcile  our  model  with  other  model.';  which  are  esiertially 
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spatial  in  outlook  (e.g.,  the  holographic  model  of  cerebral 
encoding  proposed  by  Pribram  (1971)  although  we  are  iware 
that  quasi-coherent  alpha  activity  could  provide  th.  phase 
coherent  reference  which  is  often  required  by  such  models. 

We  would  place  great  emphasis  upon  the  need  for  making  a 
careful  distinction  between  the  use  of  the  term  "phiise" 
in  referring  to  the  encoding  of  the  sensory  data  of  visual 
space  (e.g.,  Pollen,  Lee  and  Taylor,  1971)  and  our  use  of 
"phase"  as  a parameter  of  the  time  domain.  Having  cilready 
utilized  the  time  domain  to  account  for  the  encoding  of 
spatial  phase  information  (Pollen  e^  aj^.  ibid . ) , Pollen  and 
Trachtenberg  (1972)  are  led  to  regard  alpha  as  noise.  They 
comment  that  "even  if  it  were  assumed  that  all  cells  in  the 
visual  cortex  were  rhy thn^ically  excited  by  alpha  activ’ity 
and  in  phase  with  each  other,  simple  cell  responses  to  visual 
stimuli  in  different  parts  of  the  receptive  field  would  ar- 
rive at  the  complex  cells  with  diffeiert  latencies  and  thus; 
at  different  periods  of  the  alpha  cycle.  The  alpha  activity 
would  facilitate  some  inputs  and  inliibit  others  in  a non- 
predictable  way--a  noise  process.  Thus  alpha  block, whether 
by  a central  active  inhibitory  mechanism  or  by  a more  gen- 
eral desynchronization  of  rhythmic  activity,  might  serve  to 
reduce  a neural  noise  level."  In  contrast,  we  believe  that 
our  alpha  phase  modulation  model  provides  a viable  alternative 
explanation  for  the  alpha  blockage  observed  during  the  giving 
of  visual  attention  to  ncnuniform  surfaces.  According  to 
our  phase  modulatior;  mode:  1 the  d«“synchron  i za  t i or  c.f  the 
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coi'tical  response  under  these  circumstances  may  represent 
the*  complex  phase  modulation  of  a large  cortical  area  in 
response  to  complex  spatial  and  temporal  patterns  of 
sensory  input.  These  diverse  phase  adjustm.ents  (re-settings?) 
of  the  cellular  "alpha"  activities  would  greatly  increase 
the  phase  variance  of  the  cortical  cell  pool  and  thereby 
reduce  the  magnitude  of  tlie  potentials  under  a scalp  elec- 
trode. The  greater  alpha  amplitude  in  the  absence  of 
patterned  visual  contact  is  explained  by  the  greater  phase 
coherence  of  the  cortical  cell  pool  during  simple  following 
of  the  quasi-coherent  thalamic  pacemaker  activities. 

SUMMARY 

Although  simple  time  averaging  of  evoked  potentials 
confounds  and  conceals  the  influence  of  alpha  phase  contin- 
gency in  the  generation  of  the  AVER,  it  is  possible  to 
analyze  the  role  of  phase  contingency  through  the  use  of 
alpha  phase  contingent  phase  probability  measurements.  The 
resulting  analysis  is  interpreted  as  supporting  an  alpha 
phase  modulation  hypothesis  of  sensory  encoding/decoding 
in  the  cortex.  According  to  this  hypothesis,  the  sensory- 
stimulus  acts  to  phase-shift  the  activity  of  the  cortical 
receptors  away  from  the  simple  follov/ing  cf  thalamic  pace- 
makers. The  discrepancy  between  cortex  and  pacem.aker  is 
regarded  as  the  basi'-.  of  ircrebral  e*ncodirui  of  sensory  sig- 
nals. 
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A NONLINEAR  MODEL  OF  THE  HUMAN  ELECTROENCEPHALOGRAPHIC  SIG- 


NAL  DURING  PERIODIC  PHO^T I C STIMULATION 

i 

This  portion  of  the  Biocybernetics  Project  is  concerned 
with  the  development  and  testing  of  a mathematical  model 
that  is  effective  in  predicting  the  F,EG  changes  result  ina 
from  periodic  photic  stimulation. 

There  is  little  ne«’d  to  justify  the  development  of  a 
mathematical  model  for  cliaracterizing  a complex  set  of  empiri- 
cal behavior.  A good  model  is  an  essential  part  of  an  ef- 
fective computerized  bioc;ybernetic  schem.e  for  the  prediction 
and  control  of  actual  behavior.  If  the  human  brain  poten- 
tials are  to  be  used  in  a biocybernetic  system  as  estimators 
of  brain  states,  it  is  important  to  discover  the  simplest 
adequate  model  for  characterizing  the  phenomena  of  interest. 
Inefficiency  in  the  model  will  limit  the  amount  of  control 
achievable . 

The  human  electroencephalogram  (EEG)  exhibits  various 
complex  and  poor ly-under.stood  responses  to  repet itiv'e  flash 
stimuli  delivered  to  the  eye.  In  the  present  investigation, 
these  EEG  changes  are  characterized  in  terms  of  the  simplest 
nonlinear  input-output  model  having  similar  properties  of 
response.  Quantitative  testing  and  evaluation  of  the  charac- 
terization are  achieved  through  the  use  of  dinital  simulations 
of  the  model's  behavior  and  the  use  of  digital  analyses  of 
experimental  EEG  data.  "’he  fit  lietween  the  m.ridi  1 ' s predictions 
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and  the  EEG  is  good,  particularly  for  periodic  stimulx,  and 
extends  over  a relatively  wide  range  of  frequencies.  The 
model  is  reasonably  consistent  with  what  is  known  about  the 
physiology  of  alpha  wave  generation  and  it  is  able  to  unify 
several  apparently  disparate  phenomena  (i.e.,  they  are  seen 
as  derivable  from  the  sarnie  nonlinear  sy’stem)  . The  model  is 
also  consistent  with  the  sampled-data  model  (Vossius,  19fil; 
Young  and  Stark,  1962;  etc.)  for  relationships  between  sac- 
cadic eye-movements,  EEG  activity,  cortical  excitability , 
and  visual  perception. 

As  preparation  for  the  development  of  a model  of  EEG 
alpha  generation,  a large  body  of  relevant  literature  is 
examined,  including  work  on  the  EEG  as  an  autonomous  sig- 
nal, the  EEG  response  to  periodic  photic  stimulation  (flicker), 
the  time-averaged  response  to  single  stim.uli  (average  visual 
evoked  potentials),  thalamo-cortical  neurcj-liysiology , thalamic 
and  cortical  excitability  cycles,  visual  reaction  time,  sac- 
cadic timing,  and  short-term  visual  meniory . Of  critical 
imj'Ortance  to  this  model  is  the  observation  that  the  mechanism 
of  alpha  rhythm  generation  is  nonlinear . This  and  several 
other  key  phenomena  are  noted. 

To  characterize  these  phenomena,  a simple  second-order 
nonlinear  differential  equation  exhibiting  similar  behavior 
is  selected.  This  equation,  known  as  the  van  de  Pol  oscil- 
lator, is  used  to  model  the  EEG  response  to  periodic  visual 
stimulation.  Tn  this  model,  the  autonomous  limit  cycle 
oscillation  represents  the  alpha  rhythm;,  and  a croupled 

7'i 


1 


I 

1 

^ excitation  function  reprei^ents  visual  field  intensity,  P 

f 

perturbation  technique  foj'  analyzing  the  response  to  irri- 
puJ se  trains  is  developed,  and  numerical  verifications  are 
obtained  through  simulations  in  a digital  computer. 

EEG  recordings  from  subjects  exposed  to  sine-modulated 
or  to  stroboscopic  flashes,  on  the  one  hand,  and  simulated 
EEG  generated  by  the  model,  on  the  other  hand,  are  processed 
identically  and  the  results  are  compared . The  behavior  of 
the  nonlinear  model  matches  the  actual  EEG  responses  quite 
well  over  a relatively  wide  range  of  stimulus  frequencies 
and  experimental  parameters.  Given  the  low  order  of  the 
van  der  Pol  oscillator,  the  number  of  EEG  phenomena  which 
it  predicts  is  rather  surjjrising  (and  gratifying)  . How- 
ever, the  fact  that  the  model  does  not  predict  transient 
effects  as  well  as  it  doe.s  steady-state  behavior  suggests 
that  a higher-order  model,  such  as  a network  of  oscillators, 
is  worth  investigating. 

It  must  be  recognized  that  an  input-output  model  of 
the'  sort  employed  here  is  not  a physiological  model  because 
there  is  no  plausible  correspondence  between  the  model  com- 
ponents and  particular  physiological  mechanisms.  Much  of 
the  success  of  the  oscillator-model  used  in  the  pr<;sent 
study  is  due  to  the  fact  that  it  is  treated  as  a canon ica 1 
model,  a tool  for  classifying  phenomena  in.  their  simplest 
form.  As  such,  it  provides  a quantitative  template  for  any 
model  attempting  to  account  for  t)ie  same  pfienomena  and  it 
provides  a useful  set;  of  behavioral  classes  against  which 

so 


onti  may  check  the  performance  of  an  alternative  model. 

Until  it  has  been  proved  otherwise,  we  can  entertain  the 
hypothesis  that  the  van  der  Pol  type  of  oscillator  is 
the  simplest  nonlinear  model  that  exhibits  the  same  classes 
of  behavior  as  the  actual  process.  In  any  event,  the  model 
provides  a useful  step  in  the  direction  of  developing  a 
more  complex  and  more  capacious  model. 

This  nonlinear  model  of  the  human  EEC  signal  is  des- 
cribed in  considerable  detail  in  a Ph.D.  dissertation  by 
John  R.  Nickolls  (1977),  a copy  of  which  is  subm.itted  with 
this  report. 
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APPENDIX  I 


A Classification  of  Recursive  Modeling  MetKods 

I.  Introduction 

Closely-coupled  man-machine  systems  require  efficient  estimation  and 
prediction  of  central  nervous  system  (CNS)  activities,  such  as  eye  movements  and 
EEG  signals.  Although  the  systems  responsible  for  these  CNS  activities  are  complex 
and  time- varying,  linear  systems  modeling  techniques  can  be  successfully  used  to 
predict  CNS  states.  Model  parameter  estimates  can  be  updated  in  time,  yielding  a 
time-varying  linear  model.  Nonstationarltles  in  observed  CNS  activities  can  be 
modeled  with  time-varying  model  parameters,  or  with  a time-invariant  linear 
model  with  suitable  initial  conditions. 

The  flexibility  of  linear  systems  modeling  has  resulted  in  the  rapid  development 
of  methods  and  applications  in  many  areas,  particularly  seismic  data  processing  and 
speech  modeling.  Due  to  the  diversity  of  these  developments,  there  exists  a 
plethora  of  methods  for  estimating  the  parameters  of  linear  models  given 
input-output  data,  transfer  functions,  or  covariance  functions.  Here  we  present  a 
systematic  classification  of  exact  least-squares  modeling  procedures  that  are 
recursive  (in  model-order)  and  optimal  in  some  sense  [MKLNV].  Methods  which 
are  suboptimal  or  approximate  will  only  be  briefly  indicated.  Within  this 
framework,  we  shall  point  out  some  new  algorithms  that  have  many 
computational  advantages  over  existing  ones.  Since  we  consider  state-space, 
autoregressive  moving  - average  models,  and  the  related  ladder  realizations,  we 
shall  distinguish  the  following  three  classes  of  algorithms:  Rlccatl  or  square-root 
type  methods,  recently  developed  fast  algorithms,  and  their  ladder  forms.  While 
the  first  class  typically  requires  computations  of  0(n’^)  or  CKnh  with  n equal  to  the 
number  of  model  parameters,  the  fast  forms  only  require  operations  and  storage  of 
0(n).  Because  efficiency  is  critical  in  real-time  estimation,  prediction,  and  control, 
these  fast  algorithms  are  of  particular  interest  here. 

In  Section  II  we  introduce  the  modeling  problem  by  reviewing  external  and 
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internal  (linear)  models,  and  consider  the  types  of  observed  data  used  to  estimate 
the  model  parameters.  We  then  Introduce  a systematic  classification  in  tableau 
form  of  the  various  methods  to  be  discussed.  It  should  be  stressed  here  that  these 


least-squares  modeling  methods  can  In  general  only  determine  the  unique 
innovations  representation  (IR)  model  [K-S74].  This  is  usually  a 
finite-dimensional  linear  model  driven  by  white  noise,  conceptually  the 
innovations  or  one-step  prediction  errors.  The  output  of  the  model  represents  the 
process  of  Interest,  e.g.,  eye-movement  signals  or  EEG  signals.  There  exist  many 
such  models,  but  the  IR  model  is  the  only  one  which  is  both  stable  and  stably 
invertible.  The  inverse  process  yields  the  innovations  when  driven  by  the 
observed  data.  This  means  that  even  though  the  driving  Inputs  to  the  IR  model  are 
usually  not  known  (e.g.,  neurological  control  signals),  they  can  be  estimated  from 
the  current  model  parameter  estimates.  Thus  the  parameters  of  the  IR  model  are 
chosen  to  produce  behavior  statistically  equivalent  to  the  observed  data. 

In  Section  III  we  consider  batch  methods  that  are  best  suited  for  cases  where  the 
observed  data  is  accessed  in  blocks.  This  typically  occurs  when  the  data  is  in  the 
form  of  a covariance  function  or  system  transfer  function.  A traditional  block 
method  based  on  the  Yule- Walker  equation  [Par]  has  been  used  to  predict  the  EEG 
alpha  wave,  as  indicated  elsewhere  in  this  report.  Both  autoregressive  (AR)  and 
autoregressive  - moving  average  (ARMA)  models  are  treated,  and  the  fast  versions 
which  utilize  certain  matrix  properties  are  indicated. 

In  Section  IV  recursive  (in  time)  modeling  methods  are  discussed.  These 
procedures  access  the  Input-o'  sequentially,  and  are  known  as  on-line 

methods  in  the  control  context  [..i.J,  ^MKL].  They  are  Ideally  suited  for  real-time 
blocybernetlc  applications.  Again,  the  fast  versions  of  these  algorithms  are  pointed 
out. 

In  section  V the  new  ladder  (or  lattice)  - type  realizations  of  the  fast  algorithms 
discussed  in  Section  III  and  IV  are  introduced  [IS],  [Mo].  These  new  methods  have 
nice  stability  properties  and  good  numerical  behavior.  They  also  have  lowest 
computational  complexity  and  minimal  storage  requirements. 
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In  Appendix  A we  show  how  to  obtain  partial  realizations  of  the  Joint  impulse 
response  - matching  and  covariance  - matching  type,  that  are  both  stable  and 
minimal.  In  Appendix  B,  we  present  an  example  of  our  new  exact  least-squares 
recursions  for  ladder-form  realizations,  called  the  "pre- windowing"  case  [MDKV], 
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II.  The  Modeling  Problem 

The  many  different  types  of  linear  models  can  be  classified  Into  "external"  or 
Input/output  descriptions,  and  "internal"  or  state-space  descriptions.  We  will  first 
consider  "external"  descriptions,  which  are  sometimes  referred  to  as  transfer 
function  type  models. 

Let  us  consider  a finite-dimensional  linear  system  (FDLS)  with  Inputs  | «(•) } and 
outputs  { y(’) }.  The  outputs  represent  sampled  data,  such  as  speech  where  y's  are 
scalars,  or  seismic  signals  from  a geophone  array  where  y's  are  r-vectors.  The 
Input-output  relationship  can  be  described  by  an  autoregressive  - moving  average 
(ARMA)  model 

(2.1a) 
(2.1b) 
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where  | «)(•) } Is  a moving  average  of  a white  noise  sequence  j u(-) ) and  the  values 

{y  j y_^  ) and  { u.].  . . . , u_^  ) are  initial  conditions.  The  modeling 

problem  here  is  to  determine  the  model  parameters  a-  and  . Applying  the 
z-transform 
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to  equation  (2.1)  yields 
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(2.3c) 
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With  zero  initial  conditions  and  scalar  processes,  the  ratio  of  y(z)  and  u(z)  gives  the 
transfer  function  T(z)  - b(z)la(z)  . When  Hz)  • z^  , ki  0 then  [ w(-) ) Is  a white 
noise  sequence  and  { >(•) ) is  called  an  autoregressive  (AR)  processi  the  model  is 
referred  to  as  all-pole.  Alternatively,  when  a(z)  - z",  { y-)  ) Is  a moving  average 
(MA)  process  and  an  all-zero  model  Is  obtained. 


Turning  to  "Internal"  models  of  the  FDLS,  we  can  consider  j y{')  } as  being 
generated  from  { u(<)  } with  a suitable  initial  condition  (Xg)  by  the  state-space  model 


>,•  - H . no  (2.4) 

This  model  can  be  chosen  to  have  the  given  transfer  function 

r(z)  . H(rI-*)-'rr  - ^ . (2.S) 

Note  that  for  convenience  we  have  used  a model  driven  by  rather  than  the 
more  commonly  used  model  driven  by  [MKD]  since  they  can  be  related  [Mo], 
Given  the  transfer  function,  a simple  way  to  choose  the  matrices  { H,  ♦,  F } is 
the  "observer  canonical"  form 

♦ . Z - H , H . F - 0^]'  . (2.6) 

where  ^ denotes  transpose,  ^ [a j,  . . a„]^,  ■ [f>0-  ■ . ^ Is  the  "delay 

matrix":  Z.  j • if  (j~i  - I)  then  1 else  0 , and  «j  ^ [I,  0.  ....  0]^  is  the  first  unit 
vector.  The  state-space  model  provides  a convenient  way  of  computing  the 
covariance  function  of  the  output  process.  Even  though  the  underlying  model 
( H,  ♦.  F } or  j ) is  time-invariant,  the  output  process  j y(’)  ) Is  In  general 

not  stationary  due  to  "transients"  caused  by  the  initial  conditions.  However,  If  ♦ Is 
a stability  matrix  where  all  eigenvalues  have  magnitude  less  than  one,  then  as  i-><» 
the  transients  eventually  die  out  and  the  process  becomes  stationary.  In  the 
stationary  case,  the  covariance  is  a function  of  |i-j|  given  by 

Ryii.j)  “ H n , where  Tl  is  the  stale  covariance  matrix  as  l-*w  (see 
[DKM]). 

ARMA  models  and  state-space  descriptions  are  Just  two  different  methods  of 
representing  the  input-output  relations  of  a FDLS,  and  they  can  be  closely  related  to 
each  other  using  matrix  fraction  descriptions  (MFD's)  [DKM].  A lesser-known  cla.ss 
of  FDLS  representations  are  the  ladder  realizations,  which  are  discussed  In  section 
V and  are  also  related  to  the  ARMA  and  state-space  models. 


In  modeling  a process  as  a FOLS  driven  by  white  noise,  the  observed  data  Is 
usually  available  in  one  or  more  of  the  following  forms: 

a)  input-output  pairs  j u(*),  >(•)  ): 

b)  Impulse  response  or  related  sequences  such  as  moments,  (or  moment  estimates, 
e.g.  obtained  from  Input-output  pairs); 

c)  covariance  functions  or  second  order  knowledge  of  Inputs  and/or  outputs,  (or 
estimates,  e.g.  from  the  impulse  response). 

Batch  methods  are  used  when  data  is  accessed  in  blocks,  as  in  b)  and  c);  efficient 
methods  for  determining  model  parameters  are  recursive  In  order.  Recursive  (In 
time)  methods  are  appropriate  when  data  Is  accessed  sequentially,  as  in  a).  Model 
parameters  can  then  be  estimated  recursively  both  In  order  and  time.  Table  1 
Illustrates  the  modeling  methods  that  we  will  discuss;  they  are  divided  first  Into 
batch  and  recursive  groups,  then  by  model  class;  Riccatl  or  square-root  methods, 
"fast"  algorithms,  and  their  ladder  forms.  Within  each  class  the  name  or  code  for 
each  method  appears  along  with  some  pertinent  references. 


91 


Table  I:  Classification,  of  Least-Squares  Modeling  Methods 
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MA-part,  Cholesky.PE.SQr  [GK],Sec.lV  Fast  Cholesky.RMHC:  Sec. IV, [Mo]  | Rec.Fccdb. l adder:  Sec. IV, [Mo] 

Hoc.  Spectral  Fact.  J 


III.  Batch  Methods 


When  observed  data  is  available  in  blocks,  batch  modeling  methods  are 
convenient.  We  will  first  consider  AR  models  because  of  their  widespread  use 
[Makh],  The  recently  developed  "linear  predictive  coding"  (LPC)  speech 
compression  schemes,  for  example,  are  a direct  application  of  least-squares  fitting 
of  AR  models  [AH],  [IS],  [Wak]  (for  a survey  see  [MG]).  AR  models  have  also  been 
very  useful  in  statistics  [Par],  [BJ],  spectral  analysis  [UB],  [Aka],  and  multichannel 
geophysical  applications  [Rob],  [WR]. 

Normal  Equations 

It  is  well  known  in  least-squares  problems  that  the  parameters  of  an  AR  model 
satisfy  a set  of  linear  equations  called  the  normal  equations  (see  [K-S74],  [MG])  i 

- [1,-0,  , . . . , -a„]  - [«(.0  , . . . , 0]  . (3.1) 

An  alternate  form  is  the  Yule-Walker  equation  [Par] : 

'^fi-1  “ - ‘*n^'^n-l  " f'"! 

In  both  forms  R is  a covariance  matrix  and  the  a,'s  are  the  "predictor"  or  AR  model 
parameters.  Rf  is  the  "prediction  error"  or  innovations  covariance,  a 
non-increasing  function  of  n (typically  the  model  order).  In  speech  processing,  two 
popular  methods  of  obtaining  the  normal  equation  are  the  "autocorrelation"  method 
and  the  "covariance"  method  [MG],  but  there  exist  many  ways  of  estimating  the 
covariance  R^  [BJ],  [MDKV],  [Di].  General  methods  for  solving  such  linear 
equations  Include  Gaussian  elimination  (GE),  Cholesky  decomposition.  Householder 
transforms  [Hou],  [GGMS];  however,  they  all  require  computations  of  O(n^). 

The  Levinson-Wiggins-Robinson  (LWR)  Algorithm 
An  algorithm  that  requires  only  O(ri^)  computations  for  the  recursive  solution  of 
normal  equations  with  Toplitz  covariance  matrices  (corresponding  to  an 
assumption  of  stationarily  of  the  process)  was  first  described  by  Levinson  [Lev] 
and  later  extended  by  Wiggins  and  Robinson  [WR].  By  making  use  of  certain 


shift-invariance  properties  of  Toplltz  matrices  (the  /,J-th  entries  are  only  a 
function  of  i-j),  this  algorithm  solves  the  normal  equation  via  a set  of  recursions 
that  update  the  AR  parameters  or  the  predictor  parameters  in  increasing  order 
[Rob],  [K-S74].  The  Levinson  recursions  are  also  closely  related  to  the  orthogonal 
Szego  polynomials  [Sze],  [GS],  [KVM].  Levinson's  algorithm  plays  a central  role 
because  It  can  be  generalized  to  handle  multi-channel  data  [WR], 
multi-dimensional  or  image  processing  problems  [LKM],  nonstationary  processes 
with  "shift-low-rank"  [Mo],  [FMKL],  ladder  realizations  [MV],  [MVK]  and  ARMA 
or  state-space  models  [MKD],  [MKL],  [DKM]. 

ARMA  Models 

In  state-space  terms,  the  problem  is  to  find  a triple  {H,t,r}  such  that 
T-  - H4>*r,  where  { Tj  ) is  a given  set  of  "first  order"  data  characterizing  the 
impulse  response  of  a linear  system.  This  is  the  partial  realization  problem  [KFA], 
[DMK],  The  central  role  in  this  realization  theory  is  played  by  the  Hankel  matrix 
with  entries  H.  . > T.  . , . The  columns  or  rows  of  this  matrix  are  known  to  span 
the  state-space,  so  any  method  for  finding  a basis  is  a viable  realization  method.  Of 
particular  Interest  are  methods  for  finding  the  smallest  basis  resulting  in  minimal 
order  n realizations  [HK],  [Si],  [YT];  they  all  require  0(n'^)  operations. 

From  a transfer  function  point-of-view,  the  partial  minimal  realization  problem 
is  that  of  finding  two  relatively  prime  (matrix)  polynomials  a(z)  and  6(z)  such  that 
the  given  power  series  T(z)  matches  say  k terms  of  the  expansion  of  b(z)la{z) . This  is 
the  classical  Fade  approximation  problem,  which  in  the  scalar  case  yields 
T(z)  a(z)  - Mr)  + ( iflrrn*  in  z'\  i > k-n  } Equating  coefficients  of  z'^,  0 s i s n , 
we  get 

"nt" (3.3) 

where  is  a Toplltz  matrix  containing  the  reversed  column  ordered  Hankel 
matrix  . Note  the  similarities  here  to  the  normal  equation  (3.1)  and  to  Prony's 
method  [MG]. 
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Again,  standard  methods  could  be  used  to  solve  for  but  they  all  require 
computations  of  0{n^).  However,  If  one  takes  advantage  of  the  structure  of  the 
Hankel  or  Toplitz  matrices,  fast  algorithms  can  be  found.  Such  algorithms  have 
been  developed  (In  a coding  theory  context)  by  Berlekamp  [Be]  and  for  minimal 
realization  by  Massey  [Mass].  Multivariable  versions  have  also  been  developed 
[Mo],  [DMK].  These  recursions  are  strikingly  similar  to  Levinson's  recursions;  the 
Berlekamp-Massey  algorithm  can  also  be  related  to  orthogonal  Lanczos  polynomials 
[Lane],  [Mo].  An  alternative  method  for  obtaining  stable  partial  minimal 
realizations  is  discussed  in  Appendix  A.  It  can  also  be  derived  by  considering  a 
Gram-Schmldt  (GS)  orthogonalization  on  the  columns  of  the  Hankel  matrix  or 
the  Toplitz  matrix  [Mo],  or  more  generally  by  using  projection  methods 
[KKM]. 


Spectral  Factorization  and  Innovations  Representations 
The  problem  of  obtaining  a model  of  a process  ( y{’)  ),  given  Its  covariance 
function  or  second  order  Information,  Is  called  stochastic  realization.  We  are 
interested  In  representing  | ;y(-) ) as  the  output  of  a linear  model  driven  by  white 
noise.  In  general,  there  exist  many  such  models,  however  the  only  stable  and 
stably  invertible  model  is  the  (unique)  innovations  representation  (IR)  [K-S74]. 
The  Inverse  model  is  the  whitening  filter  that  produces  a white  noise  process,  the 
Innovations  { c(>)  ),  when  driven  by  the  observed  data.  In  discrete-time  or  time 
series  analysis,  the  innovations  are  the  one-step  prediction  errors  of  the 
observations.  If  the  process  j ji(*) } is  stationary,  the  problem  of  obtaining  the  IR 
essentially  reduces  to  one  of  spectral  factorization. 

An  efficient  method  for  obtaining  the  spectral  factors  of  S^(z) 

5y(z)  - Hr)  b{-z)  I a(r)  a{-r)  , (3.4) 

Is  given  as  a two-step  procedure  in  [DKM],  [MKD],  [Mo].  In  the  stationary  and 
scalar  process  case  considered  here,  the  truncated  or  one-sided  power  spectrum  5^(r) 
of  { y(-)  ) Is  formed  from  the  covariance  sequence.  A minimal  realization  algorithm 


Is  then  used  to  approximate  S^{z)  by  q{z)la(z) , where  lla{z)  is  the  AR-part  oi  the 
desired  model.  From 

S^{z)  - a(z)  5j.(z)  a(z"')  - o(z)?<z*’)  + g{z)a(z'^) 

- Kz)<>(r‘')  (3.6) 

it  can  be  seen  that  the  spectral  factorization  problem  for  { y{-) ) is  now  reduced  to 
the  simpler  factorization  problem  for  a MA-process  ( w(*)  } , where  the  factor  Kz)  is 
the  MA-part  of  the  desired  model.  It  can  be  obtained  by  CholesKy  and  other 
factorization  procedures.  Thus  { >(•) } is  modeled  by  the  cascade  of  the  AR  and  MA 
parts  driven  by  the  innovations,  a white  noise. 

In  the  time-domain  this  corresponds  to  a factorization  of  the  (stationary) 
covariance  (a  Tbplitz  matrix)  into  triangular  factors.  The  "fast-CholesKy" 
factorization  given  by  [Mo]  is  an  efficient  algorithm  for  stationary  and 
non-statlonary  covariance  matrices  with  "shift-low-ranK”.  It  should  be  noted  that 
many  popular  covariance  estimates  have  this  property. 
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IV.  Recursive  "in  Time"  Methods 

When  the  observed  data  Is  available  as  input-output  pairs  that  must  be  accessed 
sequentially,  recursive  modeling  methods  are  the  most  appropriate.  Many  recursive 
least-squares  methods  have  been  developed  in  the  identification  and  control  areai 
they  typically  involve  solving  Riccati-type  equations  and  have  computation  and 
storage  requirements  of  0(n’^)  and  O(n^)  respectively.  Recently  "fast"  algorithms 
have  been  developed  with  reduced  computations  and  storage  of  0(n)  using  ARMA  or 
ladder  realizations . 

An  Important  set  of  least-squares  recursive  methods  for  AR-type  models  is 
described  in  detail  in  [SLG]  and  more  recently  in  [MKL].  The  discussion  Includes 
least-squares  (LS),  weighted  least-squares  (WLS),  generalized  least-squares  (GLS), 
instrumental  variable  (IV)  and  recursive  maximum-likelihood  (RML1,2)  methods 
for  ARMA  models.  All  these  methods  solve  a Riccati  equation  that  recursively 
updates  the  inverse  of  the  matrix  appearing  in  the  normal  equation  of  the  problem. 


An  alternative  to  the  Riccati  equations  are  the  square-root  forms  discussed  for 
instance  in  [MK].  They  make  use  of  the  numerically  preferrable  orthogonal 
transformations  [Hous],  [GGMS]. 

A special  case  of  the  IV  method  is  obtained  by  using  the  n-step  delayed  outputs  as 
instrumental  variables.  This  can  be  shown  to  be  equivalent  to  a minimal  realization 
problem  given  (estimated)  covariances  R . Recall  that  in  the  given  (estimated) 
covariance  case  we  discussed  a two  step  procedure.  The  first  step  was  to  obtain  a 
minimal  realization,  or  rational  approximation  of  5,(r)  by  q{z)la{z),  or  in  the 
time-domain  of  R,  by  <3  A-'  fDKM],  [MKD],  In  matrix  notation  we  obtain 

R,  A . Q . (4.1) 


where  A and  Q are  banded  matrices  of  "band  width"  n , if  the  underlying  linear 
model  is  of  that  minimal  order.  The  first  column  of  (4.1)  corresponds  to  equation 
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where  is  the  (2,1)  block  entry  of  the  appropriately  partitioned  triangular 
Tbplltz  matrix  . The  matrix  is  the  cross-covariance  of  the  last  n observations 
and  the  same  set  n time-steps  delayed.  clearly  plays  here  the  role  of  the 
reversed  ordered  Hankel  matrix  . By  noting  that  the  Riccati-type  equation  can 
be  Interpreted  as  a recursion  for  (low  rank)  updating  of  the  Inverse  of  a matrix,  we 
see  that  the  IV  method  can  be  viewed  as  a recursive  (in  time)  updating  procedure 
for  the  minimal  realization  solution  for  a in  equation  (3.3). 


Fast  Algorithms  for  Recursive  "In  Time”  Methods 

In  [MKL]  the  development  of  "fast"  algorithms  for  the  recursive  least-squares 
methods  is  discussed  in  detail.  Efficient  recursions  for  time  and/or  order  updates 
for  AR-type  models  were  first  derived  in  [Mo].  The  basic  idea  there  was  the 
observation  that  the  matrices  encountered  in  many  least-squares  problems  have  a 
"shift  Invariance"  or  a "shift-low-rank"  property.  It  can  be  characterized  by  the 
(low)  rank  ft  of  the  shifted  difference  of  a matrix  M.  fi  [ M ~ M Z ] , where 
the  "delay"  matrix  Z was  defined  in  Section  II.  This  property  Is  generated  by  the 
fact  that  these  matrices  are  sums  of  products  of  Tbplltz  or  Hankel  matrices.  It  can 
also  be  used  to  obtain  fast  Cholesky  algorithms  for  MA  processes,  thus  obtaining 
recursive  whitening  filters  of  the  AR  type  ( e.g.  RMH5  algorithm  in  [Mo]). 

Similarly  we  can  obtain  general  LWR  recursions  In  order  and  time  for  AR 
processes,  l.e.  updating  the  MA  prediction  (whitening  filter)  parameters  a-.  A 
surprising  feature  of  the  fixed-order  recursive-ln-time  algorithms  is  that  explicit 
updating  of  the  covariance  estimate  is  unnecessary,  basically  because  the  model 
parameters  are  an  implicit  characterization  of  the  covariance.  Since  the  details  of 
these  algorithms  can  be  found  in  [MKD],[DKM],[Mo],[FMKL]  , we  shall  only  give  a 
comparison  of  the  LWR-type  algorithms,  assuming  that  the  reader  is  already 
familiar  with  the  Levinson  recursions  as  described  in  [Wle],  [WR],  or  more 
recently  [K-S74]. 

The  recursions  for  the  generalization  of  the  LRW  algorithm  for  covariance 
matrices  exhibiting  a shift  invariance  property  have  a very  similar  form  to  the 


original  Levinson  recursions.  However,  In  contrast  to  the  two  solutions  required  in 
the  LWn  algorithm,  the  so-called  forward  and  backward  predictors,  we  require  in 
general  more  solutions  for  non-Tbplltz  matrices.  It  turns  out  that  the 
"shlft-low-rank"  p of  the  covariance  matrix,  regardless  of  its  size,  is  equal  to  the 
number  of  solutions  required  in  the  recursions.  For  the  case  of  covariance 
estimates  that  can  be  written  as  products  of  two  Tbplltz  matrices  (typically 
containing  input-output  data),  the  number  of  solutions  required  is  at  most  four  in 
the  scalar  case,  and  2m+2  for  m-channel  data. 

For  combined  ARMA  models  we  can  either  attempt  to  model  first  the  AR  or  the 
MA  part  and  then  try  to  estimate  the  remaining  part  of  the  model.  In  the  batch 
methods  of  Section  III  we  discussed  ways  of  obtaining  the  AR  part  first  via  minimal 
realization  and  then  the  MA  part  via  spectral  factorization.  The  other  order  of  first 
obtaining  the  MA  part  could  be  obtained  by  working  with  (an  estimate  of)  the 
Inverse  of  the  covariance  matrix,  the  so-called  informatton  matrix  [MK]. 

The  cascaded  approach  can  be  carried  out  also  in  time  recursive  form  by 
estimating  the  AR  part  via  a (fast)  recursive  form  of  the  IV  method,  as  discussed 
above,  cascaded  by  the  fast  Cholesky  recursions  for  a MA  process  (e.g.  RMH5  in 
[Mo]).  The  only  difficulty  now  Is  that  both  parts  estimate  the  models  in  the 
so-called  controller  or  "tapped  delay  line"  realization,  a dual  form  to  the  observer 
form,  which  cannot  be  merged  by  inspection. 

The  Joint  Innovations  representation  approach  discussed  in  Section  III  and 
Appendix  A is  Ideally  suited  for  recursive  in  time  methods.  Even  though  the 
driving  Inputs  (conceptually  the  innovations)  are  usually  not  available,  they  can  be 
replaced  with  their  best  estimates  obtained  by  using  the  best  current  parameter 
estimates.  This  Is  clearly  only  possible  for  methods  with  sequential  data  access.  It 
turns  out  that  this  seemingly  "suboptimal"  approach  of  substitution  has  itself 
optimality  properties  (see  e.g.  [SLG]);  a similar  situation  occurs  in  detection  of 
unknown  signals,  and  in  the  famous  separation  result  of  linear  quadratic  control 
using  state  estimates  (KFAJ.  The  recursive  maximum  likelihood  methods  In  [SLG] 
and  [MKL]  can  be  derived  from  such  an  approach.  Once  estimates  of  current 


prediction  errors  are  obtained,  they  can  conceptually  be  treated  as  Known  data,  and 
entered  for  Instance  In  normal  equation  expressions.  The  only  problem  that  might 
arise  lies  in  theoretical  proofs  of  the  convergence  of  such  methods. 


V.  Ladder  Realizations 

In  Section  II  we  discussed  the  realization  of  a given  transfer  function  T(z)  > 
b(z)la(z)  via  transfer  function  or  state-space  type  models  such  as  ARMA,  controller, 
or  observer  canonical  forms.  If  the  roots  of  a{z)  are  known,  T(z)  can  be  represented 
by  a partial  fraction  expansion.  Using  polynomial  evaluation  formulas  we  can 
obtain  the  so  called  Jordan  canonical  or  parallel  decomposition  form  (even  for 
matrix  transfer  functions)  [Mo].  The  Jordan  form  has  the  nice  property  that 
stability  can  be  checked  by  merely  inspecting  the  magnitude  of  the  roots  . Finding 
the  roots,  however,  is  numerically  sensitive. 

The  ladder  (or  lattice)  canonical  realizations  provide  a very  promising 
alternative.  They  also  have  the  property  that  stability  and  even  minimum-phase 
can  be  checked  by  inspecting  the  PARCOR  or  reflection  coefficients  [IS],  [Wak], 
[MG],  [Mo],  [Cla2].  In  contrast  to  the  methods  for  finding  roots  of  polynomials  this 
requires  only  a finite  algorithm,  the  Schur-Cohn  test  for  stability.  This  is  actually 
equivalent  to  the  Levinson  or  orthogonal  Szego  polynomial  recursions  performed  in 
decreasing  order  on  0^  or  a-(z) , [K-S74].  Given  the  stationary  covariance  matrix  , 
l.e.  second  order  information,  the  a^’s  and  the  reflection  coefficients  can  be 
computed  via  the  LWR  algorithm.  From  a stochastic  process  point  of  view  we  can 
identify  these  coefficients  with  the  partial  correlation  (PARCOB)  coefficients  or 
singular  values.  They  also  have  physical  significance  in  the  scattering  theory  of 
waves  [IS],  [Wak],  [K-S74],  [MV],  [MVK]. 

The  Levinson  algorithm  can  actually  be  carried  out  using  only  the  reflection 

coefficients  as  parametrization,  since  the  inner  product  of  a vector  r ^ (r| r,-)^ 

and  the  coefficient  vector  S',)  can  be  obtained  as  the  current  output  k-  of  a ladder 

realization  driven  by  a previous  input  sequence  containing  (rj r^}  . Similarly 

we  can  translate  other  algorithms  for  AR  models,  such  as  the  various  generalized 
Levinson  algorithms  [Mo],  [DKM],  [LKM],  Appendix  A,  into  their  ladder  form 
equivalents  as  in  [MVK],  [MV],  Appendix  B.  These  forms  are  of  interest  by  virtue 
of  their  stability  properties  and  their  numerical  robustness  --  they  typically 
require  sample  correlation  operations.  These  forms  have  also  canonical  [Mo]  and 
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Invariance  properties  [MVK],  as  well  as  minimal  storage  requirements  for  modeling 
algorithms,  as  seen  from  a comparison  of  the  the  recursions  for  the  PARCOR  and  the 
Clj  parameters  In  Appendix  B. 


ARMA  models 

Recall  that  the  first  step  of  realizing  an  ARMA  model  in  Section  HI  was  a minimal 
realization  problem.  The  solution  to  this  MR  problem  can  be  carried  out  In  ladder 
form  by  using  a Berlekamp  Massey  (BM)  - type  algorithm.  These  recursions  are 
actually  very  similar  to  the  LWR  recursions,  as  noted  in  [Mo]i  therefore  we  can  use 
an  analogous  derivation  to  obtain  ladder  forms  for  the  BM  recursions,  as  presented 
In  [GrMo]. 

Alternatively,  the  Joint  IR  approach  explored  in  Appendix  A,  leads  (even  for 
scalar)  processes  to  the  theory  of  multichannel  ladder  realizations  of  the  AR  type 
discussed  above.  Since  we  embedded  the  underlying  ARMA  model  In  a two  channel 
AR  model,  the  IR  model  will  again  be  of  order  2n  , l.e.  non-minimal.  This  would  also 
hold  for  a ladder  realization. 

Minimal  models  can  be  obtained  by  merging  the  AR  and  MA  parts  in  the  observer 
form.  It  is  also  possible  to  obtain  a minimal  rational  ladder  form  [MG],  [Mo].  The 
basic  idea  in  state-space  terms  is  to  add  a suitable  input  matrix  (F)  or  output  matrix 
(H)  to  a ladder  form  realization  of  the  AR  part  of  the  model;  this  is  possible  since 
the  ladder  forms  are  controllable  ( or  their  dual  observable  ) [Mo]. 

The  second  step  of  the  stochastic  realization  procedure  of  Section  HI  requires  a 
spectral  factorization  for  the  determination  of  the  MA  part.  As  indicated,  we  need 
to  determine  the  triangular  factors  of  the  (banded)  covariance  matrix  of  the  MA 
process.  They  can  be  obtained  from  the  Cholesky  factors  or  the  RHS  of  the  LWR 
recursions.  Similarly  it  is  possible  to  obtain  the  ladder  realizations  of  the  MA  part, 
since  the  fast  Cholesky  recursions  "by  rows"  have  the  form  of  a state-space 
equation  with  a dynamic  matrix  4>  that  has  precisely  the  same  form  as  the  4>  matrix 
of  a feedback  ladder  form  in  state-space  notation  [Mo].  As  for  the  LWR  recursions, 
there  similarly  exists  a ladder  form  of  the  fast  Cholesky  algorithm  that  requires 
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only  reflection  coefficients  as  parametrization. 

Ladder  Forms  for  Recursive  "in  Time"  Methods 
The  ladder  forms  for  exact  least-squares  solution  to  AR  modeling  have  been 
developed  In  [MV],  In  Appendix  B,  we  shall  present  the  simplest  one  of  the  many 
possibilities  of  such  ladder  forms,  the  "prewindowing"  case  [MDKV],  It  Is 
interesting  to  note  that  the  partial  correlation  coefficients  are  computed  as  sample 
cross  correlations  between  the  "forward"  and  "backward"  prediction  errors  as 
expected  from  the  stochastic  derivatif.,s  of  the  ladder  forms  [Wak],  [Mo],  [SKM], 
The  ladder  forms  of  the  GLS  and  BMLl/2  methods  discussed  in  [SLG],  [MKL]  and 
Section  IV  can  be  obtained  by  embeddnsg  the  models  in  an  appropriately  augmented 
AR  model  as  In  Appendix  A . The  IV  ntethod  led  to  nonsymmetrlc  Rlccatl  equations, 
therefore  the  fast  versions  also  require  a nonsymmetrlc  form  of  the  LWR 
recursions.  However  It  is  clear  that  these  recursions  are  then  of  the  BM  type  since 
this  algorithm  also  works  with  nonsymmetrlc  (though  triangular)  Tbplltz 
matrices.  Therefore,  we  could  obtain  "nonsymmetrlc"  ladder  forms  of  the  type 
given  in  [GrMo].  Although  the  final  algorithms  of  these  ladder  forms  are  simple  to 
implement,  the  exposure  of  the  "shift-invariance"  is  in  general  nontrivial  [MV]. 

Our  preliminary  experience  with  the  numerical  properties  of  these  algorithms 
has  been  very  encouraging;  in  general  ladder  realizations  are  superior  to  direct 
forms  for  computing  estimates  of  the  coefficients  of  a(z)  and  IKz)  . Stochastic 
approximation  or  gradient  type  methods  using  ladder  forms  can  be  obtained  easily, 
e.g.  [SV];  however,  they  have  drawbacks  similar  to  other  stochastic  approximation 
methods,  especially  for  covariance  matrices  with  extreme  eigenvalue  distributions. 
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Conclusions 


We  have  presented  a number  of  modeling  methods  in  a framework  which 
includes  the  new  ladder-form  realizations.  Of  the  three  types  of  procedures 
considered,  e.g.  Riccatl  or  square-root  methods,  fast  methods  utilizing  matrix 
shlft-lnvarlance  properties,  and  ladder-form  methods,  the  latter  two  are  most 
suited  to  problems  requiring  efficiency.  In  particular,  the  recursive  In  time 
versions  lend  themselves  to  on-line  or  real-time  applications  because  the 
input-output  data  Is  accessed  sequentially  one  sample  at  a time.  Also,  new 
parameter  estimates  are  available  at  each  sample  time,  which  facilitates  on-line 
control  problems. 

The  ladder-form  methods  also  have  the  desirable  property  that  their  stability 
can  be  checked  merely  by  inspection  of  the  reflection  coefficients.  In  addition, 
they  are  numerically  robust  since  the  major  operations  are  sample  cross 
correlations.  They  also  have  minimal  storage  requirements  for  least-squares 
modeling  methods.  The  structure  of  the  ladder-type  realization  suggests  that  the 
reflection  coefficients  may  have  physical  significance  In  the  process  being 
modeled.  For  Instance,  in  speech  models,  the  PARCOR  coefficients  correspond  to 
acoustic  wave  reflection  coefficients  In  the  vocal  tract. 
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Appendix  A:  Stable  Partial  Minimal  Realizations 
In  this  appendix  we  show  how  to  obtain  stable  partial  minimal  realizations  of 
the  Joint  Impulse-response  and  covariance-matching  type.  It  will  turn  out  that  we 
can  obtain  an  ARMA  model  by  Imbedding  It  In  a two  channel  AR  modeling  problem. 
Given  an  ARMA  model  as  represented  by  the  difference  equation  of  section  II,  we 
can  rewrite  It  as  (let  q • n) 

y,  * “l.'»i-l  + + - *1  '<,-1  - - - V,  , (Al) 


or  B^y,  ~ = b(,u,  , where 

a’  - I 1 . a,  . . ■ ■ . <i„  1 . y/  - t >,  . ■ 

fo.  *1 - f«, . • 

Now  consider  the  following  augmented  equation 


J’l-n  ^ • 
yi-n^  ■ 


(A2) 

( is  the  first  unit  vector).  This  equation  can  be  Interpreted  as  an  AR  model  for 
the  Joint  process  { y,  u } [Mo],  since  the  RHS  Is  equal  to  the  Joint  Innovations  of 
{ y,  U } , since 
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Deterministic  Case 

We  first  consider  the  deterministic  case  where  we  are  given  Impulse  response 
data  or  the  Markov  parameters.  Writing  the  input/output  relationship  In  matrix 
notation  (see  sec.  Ill)  yields 
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where  T„  is  a lower-triangular  and  ^•p  is  a full  Toplitz  matrix  of  the  Markov 
parameters  (lip  is  the  column  reverse  ordered  Hankel  matrix).  Letting  r-*®,  we  get 
the  normal  equation 
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Stochastic  Case 

From  a stochastic  process  point  of  view  we  can  express  the  normal  equation 
associated  with  the  augmented  AR  model  as 
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We  can  solve  for  the  normal  equation  of  i 

^ J ’ I 5 (A7) 

The  equations  (A5),  (A6),  and  therefore  the  non-Toplltz  equations  (A7)  (!)  can  he 
solved  recursively  with  the  LWR  algorithm.  Note  that  If  - 0 , the  minimal 
order  n - Jk.  We  could  bring  equations  (A5),  (A6)  into  a more  familiar  form  by  the 

Interleaving  permutation  (1,3,5 2n-  I,  2,  4, 6 2n42),  cf.  [MDHV],  to  convert 

the  two-process  covariance  matrix  into  a n by  n block  Toplitz  matrix,  with  2 by  2 
blocks,  however  the  LWR  algorithm  clearly  applies  to  both  representations  with 
suitable  modifications. 

Thus  we  have  shown  that  the  Joint  impulse  response/covariance  matching 
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problem  Is  equivalent  to  solving  a set  of  normal  equations  associated  with  a two 
channel  AR  modeling  problem.  Since  the  predictor  for  the  Joint  process  Is 
triangular  and  minimum  phase,  the  denominator  of  the  underlying  ARMA 
model  Is  also  minimum  phase  and  therefore  stable,  (for  all  k). 

Equations  (A5)  and  the  elegant  stability  proof  were  actually  first  obtained  by 
Claerbout  [Clal]  via  a least-squares  rational  approximation.  The  connections 
between  the  joint  Innovations  representation,  the  augmented  normal  equations, 
and  the  Hankel  matrix  were  pointed  out  in  [Mo]  and  also  In  [MDKV],  [MKD], 
[DKM],  where  algorithms  were  given  to  solve  equations  of  the  type  seen  In  (A6) 
and  (A7).  For  the  special  case  where  H has  a "shift-low-rank"  of  one  of  the  type 
(B5),  called  the  post-windowing  method  In  [MDKV],  a Levlnson-type  algorithm 
was  given  recently  by  [MR].  The  stability  property  of  the  AR  model  was  proved 
there  using  a somewhat  more  complicated  Lyapunov  technique. 
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Appendix  B:  LS-Recursions  for  Ladder  Forms 

The  Prewindowing  Case 

Given  a series  of  observations  0<tiT],  where  ( y{-)  ) can  be  m vectors,  we 

wish  to  find  the  least-squares  one-step  predictor  of  order  p parametrized  by  the 

(matrix)  coefficients  , p)  ■ We  can  define  many  different  squared 

error  criteria  y,  for  Instance  as  a function  of  s and/  in 

/ 

^p,T  “ S ■ ‘p,<  ■ - 




An  obvious  choice  from  an  innovations  point-of-view  is  (J-0./-D,  the 
"pre-windowing"  case  [MDKV],  It  s = p and  f = T the  so-called  "covariance" 
method  is  obtained,  and  if  i ■=  0 and  / = 7"  + p we  get  the  pre-  and  post-windowed 
case  OT  the  "correlation"  method  [MG],  The  total  squared  error  can  be  expressed  as 
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Thus  the  problem  of  determining  A 7.  by  minimizing  £„  7-  leads  to 
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Although  ^p-y  is  not  Toplitz,  it  still  carries  a certain  shift-invariance  structure, 
given  by  the  following  identities 

«p7-  - flpj-.,  + y[TT-p]  yiTT-py 
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Define  the  backward  predictor  Bp  y and  the  smoothing  errors  Cpj 

B\,r  «p,r  0 , . . . 0 . ; c\  y Rp  y : y\T.T-p)  (B6) 

Then  the  forward  and  backward  prediction  errors  (innovations),  <p  y . and  y , and 
an  auxiliary  scalar  >p  y can  be  defined  by 

’’'p,T'  ’’'p.r^  ■ t -^p.r  - ®p,r  ' 


Order  Update  Recursions 


Using  the  three  shift-invariance  identities  for  (B6)  and  using  some 

symmetry  properties,  the  order  update  recursions  for  ^pX  < ^pT  ’^pT  ' ^ p T ’ 
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The  order  update  recursions  are  very  similar  to  the  multivariate  version  of  the 
Levinson  algorithm,  and  a similar  set  of  recursions  for  time-update  can  also  be 
obtained  [MDKV],[Mo], 

Ladder  Type  Realization 

Premultiplying  the  above  equations  by  ylTT-p+l]  , we  obtain  the  following 
order  update  recursions  for  7. , r„  7. , 7. 

Pi*  Pi*  Pi* 


‘p.i.r  ■ ‘p,r 


^\,T^"p,r-i  '■p.T’-i 

^'p,r^  ‘p,r‘p,r 


y..r  * r' 


p*i,T"  p*i,7  >»i,r 


The  "Kalman  gain"  A;.,  7.  is  obtained  from  [MV]  as  follows 
Pi* 

^p,/’4i  ‘ ^p,r  + '■p,r‘Vj‘i  ^ <‘"’'p-i,r^ 

and  the  reflection  or  PARCOfl  coefficients  are  obtained  by 

K*.  ,,  ^ K -r  ^ ‘ r ; K'  t ~ K'  t fi  '-T  1 

I 7 1,1  1,1  ’ 1,1  1,1  1,1-1 


(BIO) 


(BM) 


Tlie  initial  conditions  aro  tiivon  by 


for  p>T  : 
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As  the  dual  to  the  stochastic  forms  in  [IS],  [Wak],  [Mo],  [SKM],  equations 
(B8)-(B11)  are  a complete  set  of  order  and  time  update  recursions  to  obtain  the 
exact  least-squares  ladder  form  predictor,  which  is  shown  in  Figure  1. 
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1 ^ 
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Figure  1 . Ladder  realization  of  exact  onc-slcp  least-squares  predictor. 


The  recursion  (BIO)  computes  the  sample  cross-covariance  of  the  forward  and 
backward  innovations,  using  the  optimal  weighting  1/(1-Y . , .),  compared  to  other 
suboptimal  schemes  [SV],  In  the  scalar  case  y>0  if  )/q=0,  or  in  general  if  f 1 . 

since  0<V  [ML)KV],  If  m>l,  wc  require  Tifum.  These  singularities  can  be 
avoided  by  including  a priori  estimates  of  the  covariance  i',,,  or  equivalently 
including  a weighted  norm  of  the  predictor  a,  in  the  error  criteria  Several 

such  modifications  heivt*  li^*on  provon  useful  ui  actual  iiuploiueulations. 
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